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this  work  is  to  make  the  Recursive-IDPCM  more  efficient  without  increasing  the  coding 
complexity  by  adaptive  schemes. 

The  details  of  adaptive  schemes  are  discussed.  Several  algoriths  of  subimage  activity 
classification  are  proposed  and  evaluated.  Optimum  quantizers  are  designed  according 
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result  of  the  computer  simulation  demonstrates  a  high  compression  ration  with  a  good 
subjective  reconstructed  image  fidelity  without  drastically  increasing  computation  time 
overthat  required  for  Recursive-IDPCM  has  been  achieved. 
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preserve  diffraction-limited  information  in  the  resulting  image,  both  in  the  point 
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This  work  involves  further  development  of  the 
Recur sive-IDPCM  image  data  compression  method.  The  goal  of 
this  work  is  to  make  the  Recursive-IDPCM  more  efficient 
without  increasing  the  coding  complexity  by  adaptive 
schemes . 

The  details  of  adaptive  schemes  are  discussed. 
Several  algorithms  of  subimage  activity  classification  are 
proposed  and  evaluated.  Optimum  quantizers  are  designed 
according  to  data  sources  to  minimize  quantization  error. 
Difference  data  are  quantized  adaptively  based  on  the 
objective  measure  of  subimage  activity  at  each  recursion. 
The  result  of  the  computer  simulation  demonstrates  a  high 
compression  ratio  with  a  good  subjective  reconstructed 
image  fidelity  without  drastically  increasing  computation 
time  over  that  required  for  Recursive-IDPCM  has  been 


CHAPTER  1 

INTRODUCTION' 

With  the  continuing  growth  of  modern  communication 
technology,  the  demand  for  image  transmission  and  storage  is 
increasing  rapidly.  Limitations  on  bandwidth  for 

transmitting  image  data  over  digital  communication  channels 
have  led  researchers  to  look  for  efficient  and  feasible 


methods  to 

compress 

image  data 

under  a 

given  amount  of 

distortion  . 

3ef ore 

discussing 

various 

data  compression 

algorithms,  it  is  necessary  to  examine  how  image  data  can  be 
compressed . 

The  key  point  is  that  a  large  amount  of  redundancy 
exists  in  a  uniformly  sampled  image.  This  redundancy  can  be 
explored  both  in  the  frequency  and  spatial  domain. 
Removing,  or  at  least  reducing  the  redundant  information 
from  the  original  image,  is  the  goal  of  data  compression. 
Figure  1  contains  the  cosine  transform  of  a  32  by  32  image. 
It  can  be  seen  that  the  energy  in  the  transform  domain  tends 
to  be  clustered  into  a  relatively  small  number  of  transform 
samples  in  the  low-frequency  region.  To  achieve  data 
compression,  transform  samples  of  low  magnitude  can  be 
discarded,  or  coarsely  quantized  in  a  digit  a’l  transmission 


INTRODUCTION 


This  is  the  final  report  for  the  Grant  Year  1984-1985  on  Grant 
Number  AFOSR-8 1-0 17  0.  During  this  year  two  pieces  of  research  were 
concluded  under  sponsorship  of  the  Grant.  The  two  research  activities  are 
included  as  two  separate  divisions  of  this  research  report.  The  research 
activities  are  as  follows: 

1.  Adaptive  Recursive  Interpolated  DPCM  for  image  data  compression 
(ARIDPCM).  A  consistent  theme  in  the  research  supported  under  Grant  Number 
AFOSR  under  Grant  AFOSR-8 1-0 170  has  been  novel  methods  of  image  data 
compression  that  are  suitable  for  implementation  by  optical  processing. 
Initial  investigations  led  to  the  IDPCM  method  of  image  data  compression. 
Subsequent  modifications  to  the  IDPCM  algorithm  were  investigated,  and  led 
to  the  improvement  of  the  algorithm  for  both  image  quality  and  decreased 
bit  rate.  The  ARIDPCM  version  of  the  algorithm  is  the  most  powerful  form 
yet  developed.  The  ARIDPCM  algorithm  was  fully  demonstrated  in  a  M.Sc. 
thesis  by  Mr.  Eng  Yuan  Fu.  This  thesis  is  included  as  the  section  of  this 
report  entitled  "ARIDPCM." 

2.  Deblurring  images  through  turbulent  atmosphere.  A  common 
problem  in  astronomy  is  the  imaging  of  astronomical  objects  through  the 
turbulence  caused  by  microscale  fluctuations  of  the  atmosphere.  The 
microscale  fluctuations  limit  the  resolution  of  any  object  by  ground-based 
telescope,  the  phenomenon  of  stars  "twinkling"  being  the  most  commonly 
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observed  form  of  this  degradation.  This  problem  also  has  military 
significance  in  limiting  the  ground-based  observation  of  satellites  in 
earth  orbit.  As  concerns  about  SDI  arise,  the  observation  of  Soviet 
satellites  becomes  more  important,  and  this  observation  is  limited  by 
atmospheric  turbulence.  Research  has  been  conducted  under  Grant  Number 
AFOSR-81-0170  to  study  the  ability  to  use  techniques  such  as  optical  and 
digital  computation  in  removing  the  blur  from  astronomical  images  caused  by 
turbulence.  This  research  culminated  in  a  Ph.D.  thesis  by  Dr.  Karen  West. 
This  thesis  is  included  as  the  section  of  this  report  entitled  "Correcting 
Images  for  Atmospheric  Turbulence." 

TECHNOLOGY  TRANSFER 

We  are  pleased  to  report  that  the  work  sponsored  as  ARIDPCM  has  had 
immediate  transfer  to  a  major  government  activity.  The  ARIDPCM  method  has 
been  utilized  in  a  prototype  system  built  by  a  contractor  for  a  special 
study.  The  success  of  this  prototype  has  led  to  the  designation  of  ARIDPCM 
as  a  standard  format  for  a  special  system.  This  special  system,  the 
contractor,  the  sponsorship,  and  all  details  of  the  study  are  classified, 
and  Special  Access  Requirements  apply  to  the  detailed  information. 
However,  AFOSR  sponsorship  of  Grant  Number  AFOSR-81-0170  are  essential  to 
the  successful  transfer  of  this  research  into  a  problem  of  national 
significance  and  need. 
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system,  without  serious  loss  of  resolution.  In  the  spatial 
domain,  a  large  number  of  pixels  possesses  low  inter-pixel 
variance,  or  equivalently,  they  are  said  to  be  hiahlv 
correlated.  As  a  result,  the  adjacent-sample  difference  has 
a  much  smaller  variance  than  the  'original  signal.  This  fact 
is  exploited  in  redundancy-removal  for  image  data. 


1.1  Concepts  of  Data  Compression 
We  recall  the  definition  of  an  image  element  or 
pixel,  a  digital  image  gray  level  P(i,j)  which  has  been 
discretized  both  in  spatial  coordinates  and  in  intensitv. 


'Hus  the  image  may  be  considered 


as  a  matrix  whose  row  and 


column  index  is  a  point  in  the  image  and  the  c  o  r  r  e  s  po  nd  i  n  g 
matrix  element  value  identifies  the  gray  level  at  that 
point,  i.e.,  the  pixel  value. 

Two  related  concepts  necessary  to  the  understanding 
or  this  thesis  are  the  bit-rate  and  the  compression  ratio. 
A  basic  problem  in  image  data  compression  is  to  achieve  the 
minimum  possible  distortion  for  a  given  compression  rate,  or 
equivalently,  to  achieve  a  given  acceptable  level  of 
distortion  with  the  least  possible  compression  rate.  The 
compression  is  characterized  by  bit  rate,  and  it  is 
specified  as  the  number  of  bits  per  pixel,  3PP.  Then  the 
compression  ratio  is  defined  as 


^  _  Average  bit  race  of  the  original  image  data 
Average  bit  rate  of  the  compressed  data 

It  is  further  necessary  to  define  an  objective  and  a 
subjective  image  fidelity  measure. 

The  distortion  can  be  specified  either  bv  an 
objective  measure  as  the  root-mean-square  error,  RMSE,  or  by 
a  subjective  measure  such  as  the  mean  opinion  score,  which 
is  given  by  human  observers  [  1  ]  .  The  RMSE  is  used  as  an 
objective  measure  of  the  performance  of  data  compression 
methods  in  this  thesis,  and  it  is  defined  as 


RMSE 


wnere  P  .  is  the  original  pixel,  P.  .  is  the  reconstructed 
x  J  i  j 

pixel,  and  N  is  an  original  square  image  sice. 


1.2  IDPCM  and  Its  Variants 
As  mentioned  above,  many  successful  methods  of  data 
compression  have  been  developed  both  in  the  transform  and 


the  spatial 

domains  In  recent 

years 

’  2  ’  .  Among  data 

compression 

techniques  in  the 

spatial 

domain,  extensive 

efforts  have 

been  concentrated 

on  the 

Differential  Pulse 

Code  Modulation  CDPCM)  system  [  3 '  .  The  D?CM  met  hoc  is 
particularly  attractive  because  of  its  simple  design  and  its 
rapid  speed  of  operation.  which  has  made  its  use  m 


eal-time  data  compression  of 


television  signals  possible. 


I 


out  its  major  drawbacks  nave  jeer.  a  .  so  pointed  out  in  _  -  : 
(1)  it  is  sensitive  to  variances  in  image  statistics.  (  Z it 
is  highly  sensitive  to  channel  error,  and  (3)  it  becomes 


come  .ex 


i  m  o . e  m  e  n  t  when  otter  models  are  used  such  as 


two-dimensional  autoregressive  moving  average  moce.s, 


.r.teroo.atec  c  it  fere  r.  t  la  .  Pu.se  .tie 


:at3  compression  metnoa  -as  oroorset 


is  a  tamer  c e  v e  . o  p m e n  t 


the  D  ?  C  M  met  n 


A.inougr  tr. ere  are  some  simi.arities  between  :  me 
met  hoes.  IDPCM  is  different  m  two  ma’or  wavs:  .  : 


r  o  r.  -  c  3  u  s  a  .  .  an: 


r  s  i  o  r.  o  t 


image  is  obtained  by  an  interpolation  operation, 
be  introduced  in  the  following  section.  Since 
IDPCM  method  has  been  further  imoroved.  One  is 


a  m  c  i :  r  i  e  . 


IDPCM  algorithm  called  Re cu rs i v e- IDPCM  [6],  and  another 


the  A  a  a  o  1 1 v  e  Recursive— IDPCM,  which 
following  sections. 


be  exp.aim 


A  high  data  compression  ratio  and  a  low  RMS  E  have 
been  obtained  by  the  R e c u  r  s i v e -  I D P C M  method.  However, 
real  world  images  the  statistical  properties  of  image  data 
vary  from  one  to  another,  and  even  within  a  single  image, 
some  regions  may  not  be  redundant  (and  therefore  not 
compressible),  and  other  regions  may  be  very  redundant.  It 

is  desirable  to  perform  data  compression  by  perfect  I  v 


.'-•-V ’u-  W  U-.V 


matching  the  activity  of  the  data  source. 


his  natural  lv 


leads  to  the  development  of  an  adaptive  Recurs ive-ID?CM 
method.  To  apply  an  adaptive  scheme,  several  problems  have 
to  be  solved.  A  suitable  measure  of  image  activity  must  be 
obtained,  adaptive  quantization  must  be  achieved,  and  the 
image  must  be  optimally  partitioned.  The  approaches  to  solv¬ 
ing  the  above  problems  are  examined  in  in  is  thesis,  and  once 
t  e  y  are  solved,  the  goal  of  the  Adaptive  Recurs!  ve-IDPCM 
tata  compression  mecnoc  is  to  achieve  a  nigh  compression 
ratio  with  a  good  subjective  reconstructed  image  fidelity 
wit  .tout  drastically  increasing  computation  time  over  that 
recuire:  ::r  Recursive  I T  ?  C  M . 


w-  -  .  '  T, 


CHAPTER  2 


INTERPOLATED  DPCM  AND  RECURSIVE 
INTERPOLATED-DPCM 

Before  discussing  the  Adaptive  Recursive-IDPCM 
method,  it  is  necessary  to  introduce  the  IDPCM  and  the 
Recursive  IDPCM  methods. 


2.1  IDPCM 

The  IDPCM  data  compression  method  was  developed  by 
Hunt,  as  an  optical  analogy  to  DPCM.  The  IDPCM  has  some 
similarity  to  DPCM,  but  rather  than  being  based  on 
prediction  of  future  pixel  values,  it  uses  interpolation. 
The  basic  idea  of  IDPCM  is  to  interpolate  a  subsampled  image 
in  order  to  generate  a  low-frequency  version  image  of  the 
original,  and  to  quantize  the  difference  between  low-  and 
high-frequency  information  in  a  smaller  number  of  bits  than 
the  bits  required  for  quantizing  the  entire  pixel  value. 
For  example,  if  we  subsample  an  image  by  skipping  8  pixels 
in  vertical  and  horizontal  directions,  the  low-frequency 
version  image  in  Figure  3  only  requires  1.6%  of  the  bits  in 
its  original  image.  For  the  difference  image  in  Figure  4, 
its  histogram  depicted  in  Figure  5  shows  that  most  of  the 
difference  values  are  distributed  in  a  very  narrow  region 
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around  t  h e  origin.  These  difference  values  can  thus  be 
quantized  in  a  smaller  number  of  bits  to  achieve  cata 
compression  . 

In  the  IDPCM  method  [5],  the  lev-frequency  version 
of  the  original  image  is  obtained  by  the  convolution 

=  "  ?(r,s)h[(i-r),( j-s)]  13) 

J  rs 

where  h  is  the  point-spread-function  of  the  interpolator. 
The  image  P  is  constructed  from  the  subsampied  image  by 
inserting  zeroes  in  place  of  the  missing  samples  to  make  the 
image  N  x  N.  The  difference  image  between  low-  and 
high-frequency  d( i , j  )  =  ?(  i , j )  -  ? (  i , j )  is  then  quantize:. 

The  properties  of  ?  depend  on  the  interpolation  function,  h. 

The  interpolative  representation  based  on 
minimization  of  the  ensemble  mean  square  interpolation  error 
was  derived  by  Jain  [14].  However,  this  minimization  is 
only  for  stationary  data.  For  a  two-dimensional  signs.  t_ 
i  ,  j  =  0 , 1 . . . N  + 1  such  that 

i  j  ' 

r?  ?  1  ,  - ri  - T 

1  i-n , j-m  i  j  • 

The  mean  square  error  c  of  the  interpolate:  v a . u  e 

i. 

minimized  if 
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Goes  not  possess  station  an  tv 
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estimating  non-stationary  data  aoulc  be  employed  to  optimize 
the  interpolator,  but  they  might  not  be  feasible  for  a 
real-time  compression  system.  Therefore,  some  kind  of 
approximation  has  to  be  adopted.  In  the  IDPCM  method,  the 
point-spread-function  h  is  approximated  by  a  bilinear 
interpolator  kernel.  For  example,  a  7  x  7  bilinear  kernel 
is  illustrated  in  Figure  6. 

The  IDPCM  operation  is  described  as  follows; 

1.  The  original  image  is  subsampled  at  every  fourth 

line  and  every  fourth  pixel,  to  create  a  subsampled 
image.  Each  subsample  is  quantized  in  3  bits  with  a 
uniform  quantizer,  the  maximum  and  minimum 

quantization  level  being  the  maximum  and  minimum  of 
the  original  image. 

2.  Zeros  are  inserted  into  the  missing  data  values  in 
the  subsampied  image  to  give  an  image  of  the 
original  size.  This  image  is  then  convolved  with  a 
7  x  7  bilinear  interpolator  kernel  which  is 
described  in  Figure  6. 

3.  The  interpolated  image  is  subtracted  from  the  origi¬ 
nal,  and  the  differences  are  quantized  in  N  sits. 
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ine  quantization  rule  is  a  tapered  quantizer  based 
on  the  Laplacian  probability  density  :  7  ]  ,  which  was 
found  to  be  applicable  to  the  method  of  IDPCM. 

The  subsamples  and  the  quantized  differences  are 
then  used  to  reconstruct  the  image. 


The  major  attractive  features  of  the  IDPCM  method 


It  is  well  suited  for  hardware  implementation. 

It  achieves  a  high  data  compression  ratio. 

It  is  less  sensitive  to  channel  error  than  the  DPCM 
system.  In  the  DPCM  method,  because  the  prediction 
is  based  on  the  previous  values,  a  channel  error 
will  affect  not  only  the  current  prediction  accuracy 
but  also  all  future  prediction  accuracy.  However, 
in  the  IDPCM  system,  channel  error  can  only  affect 
the  accuracy  of  pixels  within  an  interpolation 
kernel  and  usually  only  one  pixel's  accuracy  because 
a  large  amount  of  the  transmitted  data  is 
differences. 

Since  the  IDPCM  is  a  non-causa  1  operation,  it  may  be 
operated  in  parallel. 


Since  the  introduction  of  this  method,  several 
developments  to  the  IDPCM  data  compression  method  have  been 
achieved,  as  discussed  below. 


Recursive  IDPCM 


A  significant  improvement  of  IDPCM,  called  Recursive 
IDPCM,  was  achieved  by  Hunt  and  C a o  [  6  ]  .  One-dimensional 
Recursive-IDPCM  is  drawn  in  Figure  7  and  discussed  as 
follows. 

In  Figure  ”  the  pixel  values  ? ( n - 2  )  and  P ( n -  2  )  are 
subsamples.  They  are  quantized  in  6  bits  and  are  used  to 
interpolate  the  middle  point  c.  Then  the  difference  between 
the  interpolated  value  c  and  the  original  value  C,  cC,  is 
quantized  to  calculate  c'c.  The  reconstructed  middle  point 
value  which  is  the  sum  of  the  quantized  difference  c'c"  and 
the  interpolated  value  Z  together  with  sub  samples  ?  f  n -  2  '  and 
?(n-2)  is  used  to  interpolate  the  pixel  values  at  points  n-1 
and  n- 1 ,  which  are  b  and  d.  The  differences  b’b  and  d’d  are 
also  quantized.  The  difference  cc’  is  defined  as  the  first 
set  of  difference,  and  the  b’b  and  d'd  are  the  second  set  of 
differences.  It  is  obvious  that  the  second  set  of 

differences  is  calculated  by  the  reconstructed  value  c'  and 
subsamples  A  and  E;  therefore,  they  will  be  smaller  than  the 
first  set  of  difference.  This  implies  that  the  required 
number  of  quantization  bits  will  be  reduced. 

In  the  two-dimensional  case  the  low-frequency 
version  of  the  original  image,  rather  than  being  obtained  by 
one  time  convolution,  is  accomplished  by  performing  the 
recursively.  The  recursive  process  is: 


interpolation 


( - ,  sub  sample  the  original  image  (e.g.,  subsample  the  image 
skipping  3  pixels  both  in  vertical  and  horizontal 
directions)  and  quantize  the  subimage  at  6  bits  '  tne 
original  image  pixel  values  are  quantized  in  S  bits,  but 
usually  human  eyes  can  only  distinguish  32-64  gray  levels,  6 
bits  are  sufficient);  (2)  insert  zeros  into  the  sub  sampled 
image  to  double  its  size  and  then  convolve  the  ODtained 
image  with  a  3  x  3  bilinear  kernel;  (3)  subtract  the 
interpolated  values  from  their  corresponding  original  pixel 
values  and  quantize  this  first  set  of  differences  (the  above 
process  is  shown  in  Figure  8);  and  (4)  the  sum  of  the 
quantized  differences  and  the  interpolated  values  together 
with  the  previous  64  bv  64  subsamples  form  a  finer  subsample 
image,  128  by  128.  The  above  process  is  repeated  until 
reaching  the  original  image  size,  512  by  512.  During  this 
process,  the  size  of  a  subsample  image  increases  from  64  x 
64,  to  128  x  128,  to  256  x  256.  By  the  recursive  scheme, 
the  subsamples  come  closer  and  closer,  and  the  differences 
between  high-  and  low-frequency  versions  will  become 
successively  smaller.  Finally  only  the  first  subsample 
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low  variant  data,  the 

second 

set  of  differences 
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igned 

no 

quant i zat ion 

bits. 

This  omission  of 

the  last  set  and  possibly  the  next  to  the  last  set  of 
difference  values  can  save  a  substantial  number  of  bits. 


In  principle,  the  above  process  can  be  repeated  many 
times  depending  on  the  subsample's  density,  and  every 
recursion  will  reduce  the  interpolated  interval  by  half. 
Therefore,  the  recursive  process  has  a  property  that  every 
recursion  can  reduce  the  bit  rate.  The  bit  assignment  of 
the  Recursive-IDPCM  between  subsamples  are  shown  in  Figure 
9.  The  overall  bit  requirement  for  a  reconstructed  image  is 

Total  bits  =  6  x  (N/m)2  +  3ND1(N/m)2  +  [  12  (N'D2  )  (  N/m  )  2  ] 

(  6 ) 

where 

N  =  original  image  size 

M  *  number  of  pixels  by  which  sub  samples  are  separ¬ 
ated 

N  ^  *  number  of  bits  for  quantizing  the  first  set 
differences 

Nqt  *  number  of  bits  for  quantizing  the  second  set 
differences 


RIDPCM  is  a  very  efficient  data  compression  method 


which  has 

achieved  a  bit 

rate 

below  C  .  4 

and  mean-square 

error  below 

0.2/®.  Also  it 

is  a 

very  simple 

algoritnm  to  be 

implemented 

for  a  real-time 

data 

compression 

system. 

However,  this  method  is  used  to  compress  an  entire 
image  without  regard  to  the  amount  of  detail  in  a n  v 

t  o 


particular  area  of  an  image. 


A  further  improvement  is 
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Figure  9.  Bit  assignment  for  one  interpolation  kernel. 

(a)  Bit  assignment  for  subsamples  skipping  3 
pixels;  and  (b)  Bit  assignment  for  subsamples 
skipping  9  pixels. 
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maks  Recursive  IDPCM  adaptive,  which  will  be  presented  in 
the  following  sections. 
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CHAPTER  3 


ADAPTIVE  RECURSI VE- IDPCM  METHOD 


In  most 

present  image 

data  t 

ransnission  systems. 

the 

sampling  rate 

is  set  on  the 

basis 

of  the 

fastest  expec 

ted 

response  from 

the  data  source 

and 

not  on 

the  basis  of 

the 

quiescent  or 

normal  value . 

The 

main 

advantage  of 

an 

adaptive  data  compression  system  is  its  ability  to  increase 
the  compression  efficiency  to  a  maximum  for  the  specified 
data  with  tolerable  loss  of  information.  The  approach  of 
the  adaptive  Recursive-IDPCM  method  is  to  divide  an  image 
into  subimages  where  a  high  bit  rate  is  required  to  quantize 
relatively  complex  subitnages  but  a  low  bit  rate  is 
sufficient  for  relatively  simple  subimages.  To  match  the 
sampling  rate  and  the  quantization  level  to  the  subimage 
data  activity  or  complexity  would  require  an  activity 
classifier.  The  statistics  of  each  class  of  subimages  is 
calculated  for  designing  proper  quantizers.  The  subimage 
classification,  the  quantizer  design,  and  the  subimage  size 
are  discussed  in  the  following  sections. 

3.1  Sub  image  Classification 
Subimages  are  classified  by  level  of  activity. 
According  to  human  visual  acuity,  three  levels  of  activity 


are  suggested  [  3 ]  :  high  detail,  low  detail,  and  average 
detail.  The  high  detail  subin ages  are  defined  as  neighbor¬ 
hoods  of  sharp  gray -lev  el  transitions  and  low -detail 
subinages  as  neighborhoods  of  smooth  gray-level  transitions. 
The  subimage  activity  is  measured  by  statistical  image 
information,  and  three  approaches  are  proposed.  Here,  we 
define  high  detail  as  class  1,  average  detail  as  class  2, 
and  low  detail  as  class  3. 


Sub image  Classification  by 
Calculating  Sample  Variance 


Generally,  simple  descriptions  of  the  waveform  are 
orovided  bv  the  Quantities  r  9  ] 


N  ,  N 

U  -  1/N  ;  X ( n ) ,  o-  -  1/(N-1)  7  ( X ( n )  -  u  ) - 


n=  i 


n=  1 


x 

(7) 


called,  respectively,  the  sample  mean  and  sample  variance. 


The  quantity  or  ,  the  square  root  of  the  sample  variance,  is 


called  the  standard  deviation.  For  two-dimensional  image 
data,  the  sample  mean  and  sample  variance  are  defined  as 


n  n 


u..*(l/nj  7  7  P .  ,  .  , 

kfit-h 


N  -  x  L  =  x 
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where  n 

1  s 

the  size  of  square  subimages 

and  ?  is  a  pi 

value. 
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Subimage  Activity  average  detail  T2  <_  r  ^  <  T 1 

L  1  o  w  detail  T  3  <_  c  ,  <  T2 


Figures  10,  11  and  12  snow  the  results  of  the  classification 
of  subimage  activity.  By  this  method,  the  quantitative 
measure  basically  agree  with  the  subjective  measure.  Most 
detailed  regions  such  as  the  girl’s  eyes,  the  f  eatr.sr ,  ar.c 
sharp  edges  were  classified  into  class  one,  and  most  flat 
regions,  like  the  background,  the  girl's  shoulder,  etc., 
were  classified  into  class  three.  3ut  this  classification 
method  has  the  following  problems.  In  some  cases  there  is  a 
wrong  classification  on  sharp  edges  between  dark  and  bright 
regions.  For  example,  if  we  assume  that  the  subimages  have 
the  type  of  patterns  in  Figure  14,  the  variances  of  these 
sub  images  are  relatively  small.  In  other  words,  these 
subimages  will  be  classified  into  class  2  or  class  3. 
Usually  this  type  of  sub  images  are  on  sharp  edges  between 


dark  and  bright  regions 

• 

Figure  15  is 
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Equation  8. 

When 

this  image 

i  s 

enlarged,  Figure  16,  we  can  see  chat  the  sharp  edge  between 
the  girl's  shoulder  and  the  dark,  background  was  badly 
reconstructed.  Thus  is  because  the  subimages  covering  this 
edge  were  not  correctly  classified.  As  a  result, 
insufficient  bit  rate  was  assigned  to  those  subimages. 
Since  the  distortion  along  sharp  edges  is  particularly 
sensitive  to  human  eyes,  the  above  classification  is  not 
satisfactory.  The  second  problem  is  that  the  calculation  of 

cannot  be  done  concurrently  during  the  recursive 
i  j 

process.  This  will  cost  much  more  machine  time  than  the 
Recursive  IDPCM;  therefore,  this  approach  is  not  desirable 
for  real-time  implementation.  The  third  problem  is  that  the 
information  of  each  sub image  activity  has  to  be  transmitted 
to  the  receiver.  This  will  cost  some  bits. 

Subimage  Classification  Based 
on  the  Variance  of  Differences 

The  subimage  activity  can  be  associated  with  the 
statistical  property  of  the  differences  between  interpolated 
and  original  values.  This  is  possible  because  the 
interpolation  error  will  be  small  when  neighboring  pixels 
have  low  variance  or  will  be  large  when  neighboring  pixels 
have  large  variance.  Now  define  the  sample  mean  and  sample 
variance  of  difference  data  as 


(9) 


uIj  =  1/(n'1}  (kdi+k,jTi 


^  i  j 


k.L  (0,2,4. . . n- 1 ) 

K  r  L  =  0 

i,j  *  (1,  n-fl,  2n+l...,N+l) 
n  subimage  size 


In  Figure  5,  Che  histogram  of  difference  data  has  a  Laplace 
distribution  and  is  around  the  origin;  therefore,  the  s-ample 
mean  is  approximately  equal  to  zero 


=  1  /  ( n-1 )  ~ 


^  di+k , j  +  1 ^ 


The  3  is  used  to  represent  the  activity  of  subimages, 


igh  detail 


T1  <  *ij 


The  activity  of  subimages  (average  detail  T2  <_  <_  T1 


low  detail 


aij  <  T2 


This  approach  has  a  satisfactory  performance  of  subimage 
classification.  In  Figure  17,  class  one  sub images,  we  can 
see  that  most  sharp  edges,  the  detailed  regions  such  as  the 
girl's  eyes  and  the  feather,  are  classified.  In  Figure  IS, 
class  three,  most  of  the  flat  regions  are  aiso  classified. 
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interpolation  is  large. 

Another  important  advantage  of  this  approach  should 
be  pointed  out.  Since  the  difference  data  contains  the 
information  of  the  subimage  activity  and  they  are  to  be 
transmitted,  it  is  not  necessary  to  transmit  the  extra  bits 
for  indexing  each  subimage  activity.  Thus,  a  slig'ht  bit 
reduction  could  be  obtained. 

Multiplication  and  square  root  operations  require 
much  longer  machine  time  than  addition  and  logic  operations. 
Unfortunately  the  above  algorithm  has  brought  a  large  amount 
of  multiplication  operations,  which  is  opposite  to  the  goal 
of  a  real-time  data  compression  system.  A  desirable  feature 
of  a  data  compression  method  is  to  have  not  only  a  high 
compression  ratio  but  also  a  fast  operation.  This  demand 
leads  to  our  third  approach  to  the  subimage  classification. 

Subimage  Classification  by  the 
Absolute  Values  of  Differences 

As  mentioned  above,  addition  and  logic  operations 
are  preferred  for  sub  image  classification.  Gimlett  [  1 0  ] 
suggested  that  "the  weighted  sum  of  the  absolute  values  of 
the  transform  coefficients,  defined  herein  as  the  activity 
index,  is  proposed  as  an  objective  measure  of  scene  business 


"for  adaptive  transform  image  coding.  This  idea  is  adopted 
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where  n  is  a  sub  image  size,  d.  .  is  difference.  T.nis 

-  J 

algorithm  performs  the  classification  also  very  well,  and 
the  results,  figures  10,  21  and  22,  also  match  our 

subjective  measure.  Compared  with  the  second  approach,  the 
machine  time  is  reduced  about  one ''third.  Hence  the  third 
aporoac.n  is  considered  a  successful  method  to  objectively 


measure  me  sue  image  activity,  and 
classify  the  activity  of  subimages. 


will  be  use; 


high  detail 


"he  activity  of  subimages  (  average  detail 

tlow  detail 


Ootimal  Quantization 


After  the  interpolation,  the  differences  will  pass 
through  a  quantizer.  The  quantized  differences  then  will  be 


transmitted  . 


A  quantizer  is  a  device  vnose  output  car. 


only  a  limited  number  of  possible  values.  Each  input  eitrer 
analog  or  digital  signal  is  forced  to  one  of  the  allowable 
output  values.  The  input  range  is  ci viced  into  a  number  :  f 


✓  / 


bins  equally  as  illustrated  in  Figure  23.  It  an  input  falls 
into  the  k  u  n  bin,  tne  outcut  is  the  value  Q,  corresoondmg 
to  the  center  of  the  Kun  bin  so  that  each  input  is  rounded 
off  to  the  center  of  the  bin  into  which  it  falls.  A  uniform 
quantizer  has  all  its  bin  widths  equal.  Nonuniform 
quantizers  allow  different  bins  to  have  different  widths. 

In  the  Adaptive  Recursive-IDPCM  data  compression 
system,  the  quantization  strategy  is  to  choose  the  quantizer 
levels  so  that  they  minimize  the  mean-square-error.  This 


error  is 

given  by 
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e„  = 

i  =  1 

X.  , 

r  1+1 

(X 

;  X 

-  Q  )2?(x)dx 

(11) 

where  x^ 

<  x?  <  ... 

<  Xn  +  1 

and  <  Q2  < 

Q3  ...  <  Qn  are 

decision 

boundaries 

and 

the  quantizer 

output  levels, 

respectively,  and  P(x)  is  a  probability  density  function  of 
differences.  Considering  Figure  21,  the  histogram  of  the 
differences,  let  F(x)  represents  the  number  of  differences 
which  nave  values  equal  to  x  .  ,  then  the  density  function 
?f  x.  )  is  equal  to 

F(X;  <  X  <  X;  <  AX') 

PCx,)  =  lim  - - - — - - -  (15; 

AX-0 

If  we  choose  AX  =  I  (because  values  of  sub  samples  are 
integers),  F(x)  is  approximately  equal  to  the  probability 


densitv 


function  ? ( x  )  . 
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optimum  quantizer  output  levels  are  determined.  Figure  25 
snows  a  halt  range  of  3  bits  quantizer  input  and  output. 

3,3  Adaptive  Recur sive-ID?CM  Process 
The  equations  of  the  Ada p t i ve-R ec u r s i ve- IDPCM  are 

Pij  =  1/4(?ij  +  n/2,  i  +  n/2  * 


p 

*  i  - 

n/2 ,  j  -  n/2 

+ 

P  . 

1  1  - 

j 

n/2,  j  +  n/2 

+  P. 
1 

interpolation  < 

- 

■> 

j  when 

1 

f 

on  interlaced 

fie 

ir« 

I 

-  • 

W  r  1  2  *  ‘ ~ 

c 

i,  j  +  n/2 

+ 

H* 

when  on  scan  lines 


Quantization 


Z(P.  .  -  P.  .) 
iJ  iJ 


reconstruction  P.  .  =  P  *  d 

1 2  i  j  i  j 


where  n  is  the  number  of  pixels  by  which  subsamples  are 
separated.  In  Adaptive  Rec ur si ve- IDPCM ,  the  subimace  3-  ce 
is  selected  equal  to  the  number  of  pixels  5  v  •  -  _ 

original  image  is  sampled.  Figure  26  snows  a  case 
there  is  a  3:1  subsampling  of  image  oixe.s  a.:r ;  — 


*  *  subsamples  to  be  transmitted 
.  =  position  of  original  pixels 
0  *  interpolated  pixels 


Figure  26.  First  set  of  interpolation  of  an  3  x  8  pixels 
subimage. 


Interlaced 

Fields 


subsamples  to  be  transmitted 
first  set  of  reconstructed  pixels 
second  set  of  interpolated  pixels 
position  of  original  pixels 


•  l gu  re 


/  . 


Second  set  of  interpolation  of  an  3  x  3  pixels 
subimage  . 


and  8:1  subsampling  of  image  lines. 


The 


st  set  ox 


line, 


interpolation  is  shown  by  arrows. 


V  i  t 

i  the  at 

o  v  e 

equations 

,  the  interpoia: 

:  i o n  ,  zhe 

sub  image  classi 

f icat ion , 

and  the 

quantization,  the 

adaptive 

data  compression  can 

be 

conducted 

in  two  schemes. 

X  •  X  n  s 

bit 

rate 

ass 

ignment  of 

each  subimage  is 

based  on 

the  measure 

of 

sub  image 

activity 

r  r4 

bits 

*  > 

first  set 

of  differences 

Ti  <  C.  . 

-  - : 

bit  rate 

4  3 

bits 

=  > 

second  set 

of  differences 

assigned  , 

r  3 

bits 

=  > 

first  set 

of  dif  f erences 

T2 < C .  . <T1 
~  ij 

to  each  ^ 

Lo 

I  " 

bits 

=  > 

second  set 

of  differences 

s  u  b  i  m  a  g  e 

r  7 

J 

Llo 

bits 

=  > 

first  set 

of  differences 

cu  <  :: 

bits 

=  > 

second  set 

of  differences 

(25) 


where  T1  <  C.  .  indicates  high  detail,  T2  <_  C.  .  <  71,  average 

1  j  1 1 

detail;  <  T2,  low  detail,  is  equal  to  a  ^  ^  (8),  or 

a.  .  (10),  or  3-  .  (12),  and  a.  .  and  3.  .  are  calculated  from 
ij  ij  i j  ij 

the  first  set  of  differences.  However  the  reconstructed 
image  does  not  snow  very  satisfactory  quality.  The  problem 
is  that  those  sub images  classified  into  class  1  or  class  2 


were  f 

ineiv  reconstructed 

( e . g  .  ,  the  feather,  the 

girl 

• 

s 

eyes, 

etc.);  however,  the 

s  u bimages  in 

the  class 

3  we 

poorly 

reconstructed  ( e  .  g  .  , 

the  hat,  girl' 

s  shoulder. 

etc. 

) . 

On  the 

other  hand,  if  we 

look  at  the 

histograms 

of  t 

he 

second  set  of  differences  in  Figures  28  and  29,  we  find  that 
a  large  number  of  difference  data  has  zero  or  nearly  zerc 
values.  This  implies  that  the  data  compression  for 
subimages  in  classes  1  and  2  is  still  inefficient.  Of 
course,  we  could  adjust  thresholds  to  classify  more 
subimages  into  class  2  and  fewer  subimages  into  class  1,  but 
the  improvement  would  not  be  very  notable.  We  noticed  that 
in  the  Recur sive-IDPCM  method,  the  differences  are 
calculated  and  quantized  at  each  recursion.  If  we  can 
decide  the  quantization  level  at  each  recursion,  the  data 
compression  will  be  efficient.  This  leads  to  the  second 
scheme  of  adaptation. 

2.  The  differences  are  quantized  adaptively  at  each 
recursion  and  bit  assignment  is  based  on  Equations  10 
and  12.  The  quantization  bit  assignment  i_5  shown  in  the 
following  relations. 


bit  rate  for  the 


first  set  of 


difference 


quantization 


3it  rate  for  the 


second  set  of 


quantization 


&  bits 


T1  <  C.  . 

—  i  i 


3  bits  (first  recursion)  T2  <  C .  .  <  T1 

-  l  j  - 


0  bits 


3  bits 


C .  .  <  T2 
1  J 


T1  <  C. 


2  bits  (second  recursion)  T2  <_  C  <  T 1 


0  bits 


C.  .  <  T2 
1 J  ~ 
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here 


assign 

Adapt! 


C.  .  can  be  ecuai  to  a.  .  (10)  or  £  .  .  (12). 

1-  *  n  1  11 

-  J  -  ~t  x  J 

When  a  su bin age  size  is  3  x  S  pixels,  possible 
meets  on  a  subinage  are  shown  in  Figure  3C . 
ve-Recursive-IDPCM  system  diagram  is  s.nown 


- -Sure 


31. 


3.1  Performance  of  the  Adaptive 
Recursive- I DPCM 

The  Adaptive  Recursive-IDPCM  data  compression  method 
has  been  simulated  on  the  PDP  il/70-IIS  image  processing 
system  in  the  Digital  Image  Analysis  Laboratory.  The  result 
shows  that  a  further  improvement  to  the  Recursive-IDPCM  has 
been  achieved.  Comparing  the  two  reconstructed  images  in 
Figures  32  and  33,  it  can  be  seen  that  in  the  reconstructed 
image  which  was  compressed  by  Recursive-IDPCM,  the  detailed 
regions  such  as  the  girl's  shoulder,  the  edge  of  the  hat 
were  jagged  and  her  eyes,  eyelash,  and  the  feather  were 
blurred;  however,  in  the  reconstructed  image  of  the  Adaptive 
Recursive  IDPCM,  this  type  of  degradation  is  much  less  in 
evidence  even  with  a  lower  bit  rate  than  used  in  the 
Recursive-IDPCM.  3ut  in  some  simple  regions  of  the 
reconstructed  image  by  the  Adaptive  Recursive-IDPCM  method, 
degradations  are  still  noticeable,  such  as  the  top  of  the 
girl's  hat.  It  is  because  of  the  not-very  obvious 
distortion  in  simple  regions  that  a  large  number  of  bits  can 


£  a 

< 


|  I 


0  3  0 


0  0 


0  0  0 


0  0  0  0  0  0  0 


3  0  4 


3  0 


0  3  0  2  0 


0  0  0  C 


0  0  0 


00  0  00000 


0  3  0  3 


0  2  0 


0  2  0  2 


0  0  0  0 


0  0  0  0 


000000000 


6  0 


0  3  0  6 


602030206 


0  0 


0  6 


0  2  0 


0  6 


0  0  0 


00000000  0 


000000000 


202020202 


000000000 


ooooooooo 


200020002 

ooooooooo 

ooooooooo 


0  0  0  0 


0  0  0  0 


2  0 


2  0 


0  2  0  2 


OOOOOOOOO 


0  2  0  2  0 


ooooooooo 


6  0  0  0 


0  0  0  6 


602020206 


Figure  30.  Sir  assignments  for  an  8  x  3  subimage. 


be  saved  tor  detailed  regions  which  are  very  sensitive  to 
human  eyes. 

The  different  partition  of  the  original  image  was 

tested  by  dividing  an  image  into  different  sizes  of 

subimages.  The  following  table  shows  the  relation  between 
the  size  of  sub  images  and  the  number  of  subimages  in  each 
class  (Table  1). 

The  smaller  the  subimage  size  is,  the  larger  the 
number  of  simple  subimages  will  be  because  the  correlation 
between  the  nearest  pixels  is  higher.  As  a  result,  the 
difference  values  will  become  smaller  and  the  bit  rate  for 
quantizing  differences  can  be  reduced.  On  the  other  hand, 
the  overall  bit  rate  was  increased  because  the  image  had  to 
be  sampled  more  densely.  However,  a  very  g'ood  reconstructed 
image  was  obtained  with  the  bit  rate  0.546  (Figure  34). 
when  the  size  of  the  subimages  is  too  large  (e.g.,  16  x  16), 

the  subject  quality  of  the  reconstructed  image  shown  in 

Figure  35  is  not  very  satisfactory,  although  a  slight  bit 
reduction  has  been  obtained.  This  is  because  there  is  not 
much  correlation  among  too  coarsely  subsampled  pixels.  For 
the  image  used  in  this  simulation,  a  3  x  3  subimage  is  a 
suitable  size.  Generally,  cne  selection  of  optimum  subimage 
size  for  adaptive  data  compression  schemes  depends  on  the 
statistics  of  the  image  data. 


One  problem  which  should  be  pointed  out  is  that 
although  the  subjective  quality  of  the  reconstructed  image 
has  been  improved  by  using  the  adaptive  scheme,  t  n  e  objec¬ 
tive  quality  represented  by  RMSE  (2)  has  not  been  achieved. 
For  example,  in  Figure  32,  Rec u r s i v e- IDPCM ,  the  3PP  is 
0.3377  and  the  RMSE  is  0.00813;  in  Figure  36,  Acaptive 
Recursive  IDPCM,  the  BPP  is  0.3561  and  the  RMSE  is  0.00907. 
For  the  same  bit  rate,  the  Adaptive  Recursive-IBPCM  has  a 
slightly  larger  RMSE  than  the  Recursi  ve-IDPCM.  But  it  has 
been  pointed  out  [11]  that  quantitative  measures  of  image 
fidelity  which  have  been  developed  are  not  perfect. 
Therefore,  our  evaluation  on  the  Adaptive  Recursi ve-IDPCM 
method  is  based  on  its  subjective  measure. 


t  * 


CHAPTER  4 


COMPRESSION  OF  COMPUTER  TOMOGRAPHIC  PROJECTION 
3Y  THE  ADAPTIVE  RECURSIVE-IDPCM 

It  has  been  pointed  out  [12]  that  the  projection 
matrix  of  computer  tomographies  contains  a  great  deal  of 
redundant  information.  Therefore,  the  projection  matrix  is 
compressible.  The  compression  of  tomographic  projections 
using  the  DPCM  method  has  been  studied  [12].  In  this 
thesis,  a  new  approach  of  compressing  tomographic  projec¬ 
tions  using  the  Adaptive  Recursive- ID PCM  is  introduced. 

A  projection  taken  along  a  set  of  parallel  rays  is 
called  a  parallel  projection,  two  examples  of  which  are 
shown  in  Figure  38,  and  a  projection  along  parallel  rays  in 
a  certain  angle  is  calculated  by  the  function 

Ps(x)  «  J  f ( x , y ) dy 

where  f(x,y)  represents  a  two-dimensional  image  pixel  value. 

A  projection  matrix  is  also  depicted  in  Figure  39 
and  will  be  compressed.  It  has  also  been  pointed  out  [12] 
that  the  amount  of  redundancy  appears  to  be  strongly 
dependent  upon  the  angle  of  projection,  the  redundancy  is 
highest  near  the  angles  of  0°,  90°  and  180°.  Because  of 
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igure  38.  Parallel  projections 


seems  to  o  e  more  sui:ao:e. 

A  102  x  IS*  projection  matrix  in  Figure  42  was 
obtained  from  a  128  by  125  image  shown  in  Figure  40,  and  it 
was  a  102  x  134  matrix  in  Figure  42.  To  simplify  program¬ 
ming,  the  scalloped  ends  of  the  projection  matrix  were 
trimmed  before  encoding,  producing  a  rectangular  matrix. 

The  obtained  rectangular  matrix  was  compressed  by  the 

Adaptive  Rec ur s i v e- IDPCM  method  discussed  in  Section  3.3. 
The  methods  of  calculating  a  projection  and  a  back 
projection  can  be  found  from  references  [12  and  13].  The 
reconstructed  image  quality  in  Figure  42  is  improved,  which 
can  be  seen  by  comparing  the  image  in  Figure  42,  with  the 

image  in  Figure  10c  from  the  reference  [12].  However,  the 

reconstructed  image  even  without  data  compression  is 
distorted,  as  shown  in  Figure  41.  As  a  result,  the 
reconstructed  image  with  data  compression  does  not  have 
satisfactory  quality  although  some  improvement  has  been 
achieved  by  using  the  Adaptive  Rec urs i v e- IDPCM  method. 


CHAPTER  5 


CONCLUSION 

In  this  thesis,  we  have  discussed  the  fundamentals 
of  data  compression  as  well  as  details  of  the  Adaptive 
Recur sive-IDPCM  data  compression  method.  In  order  to 
implement  adaptive  schemes,  several  subimage  activity 
classification  algorithms  were  tested,  and  the 
classification  using  the  absolute  value  of  difference  was 
considered  to  be  the  best.  The  optimum  quantizer  was 
designed  to  minimize  the  quantization  error  based  on  the 
nean-squar e-error  criterion.  Especially  when  the 
differences  were  quantized  adaptively  at  each  recursion,  the 
Adaptive  data  compression  method  performed  more  efficiently. 
Compared  with  the  Recursive-IDPCM,  we  have  seen  that  the 
subjective  quality  of  the  reconstructed  image  using  the 
adaptive  scheme  has  been  notably  improved.  As  is  the  case 
for  the  IDPCM  system,  the  Adaptive  Recursive-IDPCM  system  is 
also  less  sensitive  to  channel  error  than  the  DPCM  system. 
Although  the  encoding  and  decoding  complexity  are  slightly 
increased,  the  Adaptive  Recursive-IDPCM  can  still  be  easily 
implemented  for  a  real-time  system. 


It  is  mentioned  in  the  introduction  that  image  data 
compression  methods  are  basically  categorized  into  two 
classes.  One  class  is  data  compression  in  the  transform 
domain,  and  another  class  is  data  compression  in  the  spatial 
domain.  In  the  transform  domain,  many  transform  coding 
algorithms  achieve  high  performance,  small  sensitivity  to 
fluctuation  in  data  statistics,  but  their  hardware 
complexity  is  high.  In  the  spatial  domain,  the  predictive 
methods  are  generally  easy  to  implement,  but  they  are 
sensitive  to  data  statistics.  The  Adaptive  Recursive— IDPCM 
system  seems  to  nave  both  the  advantages  of  predictive 
coding  methods  and  of  transform  coding  methods.  Real-time 
implementation  of  the  Adaptive  Recur si v e-IDPCM  method  is  the 
suggested  step  for  future  research  so  that  this  method  can 


actually  be  tested. 
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ABSTRACT 


The  turbulent  atmosphere  degrades  images  of  objects  viewed 
through  it  by  introducing  random  amplitude  and  phase  errors  into  the 
optical  wavefront.  Various  methods  have  been  devised  to  obtain  true 
images  of  such  objects,  including  the  shi ft-and-add  method,  which  is 
examined  in  detail  in  this  work. 

It  is  shown  theoretically  that  shif t-and-add  processing  may 
preserve  diffraction-limited  information  in  the  resulting  image,  both 
in  the  point  source  and  extended  object  cases,  and  the  probability  of 
ghost  peaks  in  the  case  of  an  object  consisting  of  two  point  sources  is 
discussed.  Also,  a  convergence  rate  for  the  shif t-and-add  algorithm  is 
established  and  simulation  results  are  presented.  The  combination  of 
shi ft-and-add  processing  and  Wiener  filtering  is  shown  to  provide 
excellent  image  restorations. 
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CHAPTER  1 


INTRODUCTION 


Imaging  through  random  media  is  a  problem  with  applications  in 
a  variety  of  disciplines.  It  arises  in  medical  ultrasound  imaging 
where  the  random  medium  is  biological  tissue.  It  also  occurs  in  wave 
propagation  through  the  ocean  where  there  are  inhomogeneities  due  to 
temperature  variations  and  the  presence  of  biological  material.  In  tne 
atmosphere,  turbulence  arising  from  temperature  variations  creates 
random  effects  in  images  obtained  through  ground-based  astronomy. 

In  each  case,  the  resulting  images  nave  a  characterise  ic  gram;.1 
or  "speckled"  appearance.  This  presents  special  difficulties  in  image 
restoration  because  the  degradations  are  random  in  nature,  and  if  we 
think  in  the  context  of  linear  systems  theory,  we  are  limited  to 
information  about  the  average  point  spread  function  oniy.  One  method 
developed  to  counter  this  problem,  specifically  for  the  case  of 
ground-based  astronomy,  is  the  shif t-and-add  method  (Bates  and  Cady, 
1980;  Cady  and  Bates,  1980),  which  will  be  examined  in  detail  in  this 
dissertation . 

Since  the  degradations  produced  by  the  atmosphere  are  random  in 
nature.  Chapter  2  is  devoted  to  characterizing  statistically  the 
intensity  of  an  optical  wave  which  has  passed  through  the  atmosphere. 


Although  the  literature  is  by  no  means  unanimous  in  its  choice  of  an 


1 


t  £ 


E  : 


* 


► 


r 


n 


i 

%■ 

r 

* 

V 

T 

fij* 


appropriate  probability  distribution  for  intensity,  the  lognormal 
distribution  is  heavily  favored,  both  by  theorists  and  experimental¬ 
ists.  and  we  will  present  a  physical  model  and  experimental  evidence  to 
support  its  use. 

In  Chapter  3,  we  discuss  typical  modes  of  astronomical 
lrnrgmg — short  exposures,  long  exposures,  and  speckle  interferometrv-- 
and  also  review  some  of  the  major  advances  in  obtaining  more 
information  from  such  images.  Among  the  algorithms  to  be  considered 
are  speckle  holography,  Fienup's  iterative  algorithms,  Knox-Thompson . 
and  shift-and-add . 

Since  shi f t-and-add  is  such  a  simple  and  easily  implemented 
algorithm,  the  question  of  why  it  works  naturally  occurs.  This 
question  is  analyzed  in  detail  in  Chapter  w,  which  constitutes  t.ne 
major  original  contribution  of  this  dissertation.  We  will  der-.-.e  me 
combined  point  spread  function  for  atmospheric  degradation  arc 
shif t-and-add  processing,  address  the  probability  of  error  or  :: 

"gnost"  peaks,  examine  the  case  of  extended  ooject  imaging,  arc 
determine  the  rate  of  convergence  of  the  algorithm.  Our  results  :r. 
this  chapter  indicate  that  the  shi f t-3nd-add  method  applied  to  3  ser.es 
of  short  exposure  images  may  allow  diffract  ion- 1 tm: ted  information  tc 
be  preserved. 

In  order  to  test  the  performance  of  shi f t-and-add  process. no 
and  thus  verify  our  analytical  results,  we  have  simulated  atmcspner-.  c 
turbulence  degradations,  since  no  real  data  were  available.  We  discuss 
the  two  algorithms  used  for  this  purpose  in  Chapter  c.  The  first. 


J 

bevel-coed  pv  *1c  jlarerv  is  an  akonthn  which  considers  phase 

bemaiatims  onl\  .  »e  subsequent 1  v  •modified  this  scheme  to  include 
it'.  :,iie  effects  as  -ell. 

In  "maoier  o,  we  discus?  generation  of  cegraded  images  by 
:  l  .  .";  a  simulated  uncecraaed  mace  with  the  point  spread  functions 

•  c-cuted  it.  Cnapter  5.  We  then  apply  shi  ft-and-add  processing  alone 
i" :  .>  combination  of  sh;  f  t -and-add  processing  and  Wiener  filtering  to 

•  '->e  images,  obtaining  excellent  restoration  results.  As  a  caution 
acunsi  -.akmn  excessive  claims  for  the  performance  of  this  processing, 
--  mce  mat  the  simulated  images  are  noise-free;  however,  the  combined 
• ~*rr et i : a  1  results  of  Chapter  «  and  the  simulation  results  of 

-  cemcnstrate  that  me  she  ft-and-add  method  is  an  effective  and 
c,-re-f  v  easilv  implemented  accroach  to  restoring  images  corrupted  by 
v.msmenc  turbulence  and  per.naps,  other  random  degradations. 
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CHAPTER  2 


OPTICAL  PROPAGATION  THROUGH  THE  ATMOSPHERE 

We  are  all  familiar  with  the  twinkling  of  stars  as  we  observe 
them  through  the  earth's  atmosphere.  This  phenomenon  Is  due  to  random 
phase  and  amplitude  shifts  in  the  optical  path  between  the  star  and  the 
viewer’s  eye.  These  shifts  are  induced  by  atmospheric  turbulence,  a 
condition  primarily  caused  by  temperature  inhomogeneities  arising  from 
the  action  of  winds  and  from  heat  rising  from  the  earth.  Since 
refractive  index  n  is  dependent  upon  temperature,  these  temperature 
inhomogeneities  give  rise  to  regions  with  randomly  varying  refractive 
index,  which  in  turn  result  in  the  amplitude  and  phase  shifts  mentioned 
above. 

In  ground-based  astronomy,  the  viewer's  eye  is  replaced  by  a 
telescope,  and  because  of  the  finite  time  required  to  record  an  image, 
the  random  phase  shifts  produce  a  "speckle  pattern"  as  the  image  of  an 
unresolved  star  or  point  source.  (The  tern  "speckle"  is  borrowed  from 
laser  speckle  because  of  similarities  in  appearance  and  in  some 
statistical  models  of  the  two  phenomena,  and  it  refers  to  a  granular 
structure  in  the  image.) 

Since  the  amplitude  and  phase  shifts  induced  by  the  atmosphere 
are  random,  we  wish  to  characterize  them  by  their  statistical 
properties.  In  particular,  we  want  to  determine  the  distribution  of 


5 

intensity  (or  irradiance),  since  that  is  the  quantity  which  is 
typically  measured. 


2.1  Justification  of  Lognormal  Intensity 
Statistics  for  Speckle 

We  have  chosen  the  lognormal  model  for  intensity,  as  it  is  the 
one  supported  by  the  majority  of  the  literature  (Tatarski,  1961; 
Lawrence  and  Strohbehn,  1970;  Korff,  1973;  Fried,  1966;  deWolf,  1969). 
However,  this  model  is  by  no  means  unanimously  agreed  upon,  and  a  good 
bibliography  of  the  alternatives  has  been  compiled  by  Fante  (1975). 
Since  the  lognormal  distribution  is  relatively  unfamiliar,  we  have 
summarized  its  relevant  properties  in  Appendix  A.  In  this  section,  we 
will  present  a  non-rigorous  physical  argument  for  lognormal  intensity 
statistics  and  will  also  discuss  some  of  the  pertinent  experimental 


evi dence . 


2.1.1  Physical  Grounds  for  Lognormal  Intensity  Statistics 

According  to  Strohbehn  (1968),  we  will  assume  a  laminar  model 
of  the  turbulent  atmosphere  that  consists  of  a  large  number  N  of  slabs 
oriented  perpendicular  to  the  propagation  direction  and  will  derive  the 
statistics  of  the  phase  and  amplitude  of  an  optical  wave  at  the  pupil 
of  a  telescope  or  optical  system.  We  will  depict  this  model  in  Figure 
2.1,  where  s  is  the  source,  r  is  the  receiver,  is  the  amplitude  of 
the  optical  wave  at  the  source  (with  phase  assumed  to  be  zero),  A  and  \1> 
are  the  amplitude  and  phase  at  the  receiver,  L  is  the  optical 
path-length,  d^  is  the  (random)  width  of  the  division  of  the 
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optical  path,  and  is  the  refractive  index  of  the  l1^  slab.  The 
refractive  index  is  actually  time  dependent  and  may  be  written 

n  =  n  +  An . 

10  1 

where  n  is  the  mean  of  all  the  n.,  the  refractive  index  of  the 
o  1 

atmosphere  without  turbulence,  and  An^  is  the  varying  portion  of  n.  due 

to  turbulence.  Also,  several  assumptions  are  made.  First,  it  is 

assumed  that  we  are  dealing  with  line  of  sight  propagation  with  r 

located  in  the  turbulent  medium.  This  is  the  condition  of  ground-based 

astronomy  as  opposed  to  that  of  a  satellite  looking  down  through  the 

atmosphere.  Even  stronger,  we  are  considering  only  a  straight  line 

path  from  source  to  receiver;  that  is,  ignoring  contributions  from 

scatter  at  the  receiver.  The  second  assumption  is  that  the  wavelength 

is  much  shorter  than  d^  for  every  i.  This  implies  that  we  may  use 

geometric  optics.  Third,  we  assume  that  An^  is  very  small  compared  to 

n  .  More  specifically,  An./n  is  on  the  order  of  10-^.  The  fourth 
o  r  10 

assumption  is  that  turbulence  is  homogeneous  and  isotropic,  which 
allows  us  to  use  this  essentially  one-dimensional  model. 

We  now  consider  the  amplitude  A  at  the  receiver.  The  following 
derivation  is  extremely  simplistic,  and  it  is  now  thought  that  all 
distortions  of  the  optical  wavefront  (including  those  which  cause 
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amplitude  or  intensity  fluctuations)  are  due  to  random  phase  shifts. 

In  the  case  of  intensity  or  amplitude  fluctuations,  phase  shifts  which 
occur  high  in  the  atmosphere  cause  the  various  portions  of  the 
distorted  wavefront  to  travel  in  slightly  different  directions,  thus 
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8 

resulting  in  interference.  This  interference  results  in  the  observed 
amplitude  or  intensity  fluctuations  (Lawrence,  1976).  We  shall, 
however,  proceed  with  this  simplistic  view  of  amplitude  or  intensity. 

First,  we  define  R  ,  the  reflection  coefficient  for  amplitude 

i 

at  the  boundary  between  the  (i-l)C^  and  the  ith  slabs  for  normally 
incident  light: 

^  _  reflected  amplitude 

i  incident  amplitude 


ni  ~  ni  -  1 
ni  +  ni  -  1 


nQ  *  And  -  nQ  -  An£  _  , 
no  ’  Ani  +  n0  +  “ni  -  1 


Ani  -  Ani 


2no  ♦  An.  ♦  Ln.  .  j 


which  is  very  small  due  to  the  third  assumption  in  our  model,  so 
subsequently  we  shall  ignore  all  reflected  light.  Thus  if  is  the 

amplitude  of  the  wavefront  incident  upon  the  boundary,  the  amplitude  of 
the  transmitted  wave  is 


A.  *  A .  ,  -R.A.  . 
1  1-1  1  3-1 


(1  -  R.M.^ 


We  shall  denote 


M.  -  1  -  R. 
i  i 


krf  *«._  ",  i 


•v.-Kv/a'y. 
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and  recall  tha  M  is  random  due  to  dependence  on  An^  and  An  , .  So  i f 


we  start  at  the  source  s  with  an  amplitude  of  A  ,  the  light  hits  the 


first  slab,  and  the  amplitude  transmitted  to  the  second  slab  is  . 


This  process  continues  for  each  slab  so  that  the  amplitude  at  the 


receiver  is 


A  -  M  M  A 

'  N  ‘  •  ‘  lV 


Taking  the  natural  logarithm, 


logA  a  logM^  +  .  .  .  +  logM^  +  logA^ 


where  we  may,  without  loss  of  generality,  assume  A  =  1  and  drop  the 


last  term.  We  then  assume  that  the  log.T  meet  the  requirements  for 


application  of  the  Central  Limit  Theorem.  One  such  set  of  requirements 


(Papoulis,  1965),  although  not  the  most  general,  is  that  if  we  let  a 


a.  The  xi  are  independent. 


b.  /  Z  ->  00  as  N  ->  °°, 
i  =  l  xi 

°°  a 

c.  .mx,)  p.(x.)dx.  is  finite  for  some  a  >  2. 

oo  1  r  i  1  x 


Therefore,  for  N  very  large,  logA  is  (approximately)  Gaussian,  which 


implies  that  A  is  lognormal.  We  are,  however,  primarily  interested  in 


the  statistics  of  intensity  I,  since  that  is  the  quantity  typically 


measured.  Intensity  is  related  to  amplitude  by 


I  «  A* 


Vi' 


or 


logl  =  2logA 


so  roar  iogl  is  also  (approximately)  Gaussian,  and  I  is  iognorma 1 iy 
distributed. 

We  now  turn  our  attention  to  the  statistical  distribution  of 
phase  at  the  receier.  Since  the  phase  at  s  is  assumed  to  be  0.  the 
phase  'h  at  r  is  simply  the  sum  of  contributions  from  each  slab 


*  =  L 


d .  /  a. 

i  i 


where  V  is  the  wavelength  in  the  i1-  ‘  interval  which  vanes  due  to  in, 
the  fluctuating  portion  of  the  refractive  index.  We  let  a  be  the 
wavelength  for  propagation  in  the  atmosphere  (refractive  index  n&  '  in 
the  absence  of  the  turbulent  lavers.  Then  for  the  i1  laver 


n.  =  a. /A  , 
1  l  o’ 


if)  =  2 it  y  d  .  /( n  .  A  ) 
.  **,  1  10 
1  =  1 


2tt  p 

T  L  di/rV 

0  i  =  l 


But  n  *  n  +  An .  ,  so 
1  o  i 


N 


d. 

1 


?«  -L,  n„  ♦  in- 

o  1=1  0  “1 

~N  di 

n  .  1.  l  ♦  in./n„ 
o  o  1  =  i  1  o 


Fro.17!  the  third  assumption  in  our  model, 


«.  1. 


1  .  Ani  ^2 

1  ♦  in- /n  nrt  +  ^  n  ' 

10  o  o 


=  1  - 


t  2  r— —  /  d  .  ( 1  -  in  .  /  n  ) 

Vo  i=i  1  10 


-  /  d .  -  /  d  An .  • 

Vo  i-1  1  A  <  i  =  l  1  1 

0  o 


Now  we  let  k  be  the  wave  number  corresDonding  to  N  and  note  that  the 
o  °  o 

sum  of  the  is  equal  to  L,  so 
k  L  k  N 

0  Or_. 


0  o  r* 

*  =  H - T  X  dVr 

o  n'  1  =  1  1  A 


The  constant  first  term  of  this  expression  is  always  present  due  to 


propagation  through  the  distance  L  whether  or  not  one  is  dealing  with 
atmospheric  turbulence.  The  fluctuating  (image  degrading)  portion  of 
T,  which  we  will  denote  by  <t>,  is 
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We  again  assume  that  we  may  apply  the  Central  Limit  Theorem  to  the 
random  product  so  that  for  N  very  large,  i  is  approximately 

Gaussian  where 

L/n  +d) 


1  y  _  0  0 

;  =  e 


We  have  now  completed  our  simplistic  physical  justification  for 
assuming  lognormal  intensity  statistics  as  well  as  our  determination 
that  the  fluctuating  portion  of  phase  at  the  receiver  is  Gaussian,  and 
we  will  now  consider  some  of  the  relevant  experimental  results. 


2.1.2  Experimental  Confirmation  of  Lognormal  Intensity  Statistics 
There  has  also  been  considerable  effort  to  determine  the 
correct  statistical  model  for  amplitude  or  intensity  by  experimental 
measurement.  Summaries  of  these  efforts  are  available  in  Fante  (1975), 
Strohbehn  (1971),  and  Roddier  (1981).  Strohbehn  (1971)  reports  that 
one  of  the  first  careful  experiments,  a  measurement  of  the  variance  of 
log- intensity ,  was  made  in  the  Soviet  Union  in  1965  by  Gracheva  and 
Gurvich.  Their  measurements  agree  well  with  those  predicted  theoreti¬ 


cally  using  a  lognormal  model  for  small  (<  1)  values  of  0^  j .  Ochs 


and  Lawrence  (1969)  later  measured  the  variance  of  log-amplitude  and 
concluded  that  their  data  were  in  better  agreement  with  a  lognormal 
model  than  with  a  Rayleigh  model  for  amplitude,  which  would  imply  that 
intensity  is  better  described  by  lognormal  statistics  than  by  negative 
expont  'tial  statistics.  (This  Rayleigh  model  is  due  to  an  assumption 
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that  the  real  and  imaginary  parts  of  the  electric  field  are  bivariate 
Gaussian  with  zero  covariance.  See  Goodman  (1984)  for  a  derivation  of 
this  result  in  the  context  of  laser  speckle.) 

Unfortunately,  as  Heyde  (1963)  and  Barakat  (1976)  have  noted, 
the  lognormal  distribution  has  the  unusual  property  of  not  being 
uniquely  determined  by  its  moments.  (We  have  presented  a  proof  of  this 
property  in  Appendix  A.)  Therefore  as  Barakat  (1976)  notes,  it  is  not 
valid  to  predict  a  lognormal  probability  density  function  for  intensity 
or  amplitude  simply  based  on  the  measurement  of  moments  such  as 
variance.  This  casts  much  doubt  on  the  experimental  results  reported 
above  and  also  upon  many  of  the  others  recorded  in  the  literature. 

More  recently,  careful  measurements  have  been  made  of  the 
combined  telescope-atmosphere  modulation  transfer  function  (MTF)  for 
the  speckle  interferometry  process.  This  transfer  function  has  also 
been  theoretically  obtained  by  Korff  (1973)  using  the  lognormal  model 
(although  considerable  numerical  evaluation  of  his  result  is  necessary) 
and  by  Dainty  (1973)  using  a  Gaussian  model  for  complex  amplitude,  an 
approach  equivalent  to  assuming  a  negative  exponential  intensity 
distribution.  Chelli  et  al.  (1979)  have  presented  experimentally 
obtained  MTF's  for  infra-red  stellar  speckle  interferometry  and  have 
found  them  to  be  in  good  agreement  with  Korff' s  model.  Aime  et  al. 

(1979)  have  also  experimentally  determined  the  te iescope-aimosnnere  v77 
and  have  found  their  results  to  be  in  better  agreement  wr.r.  Kartr's 
model  than  with  the  Gaussian  model.  These  experiment  a.  res./  >  >.  " 
to  the  validity  of  assuming  a  lognormal  cist  r.  '  -  ••• 
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Although  the  lognormal  distribution  for  intensity  seems  to  nave 
more  theoretical  and  experimental  support,  the  negative  exponential 
distribution  is  still  widely  used  in  calculations  appearing  m  the 
literature.  This  is  due  in  part  to  the  fact  that  Goodman  (I98<1)  has 
developed  an  extensive  collection  of  analytical  results  using  this 
distribution  in  connection  with  laser  speckle.  Also,  Dainty  (198<*)  has 
noted  that  the  negative  exponential  distribution  becomes  a  better 
assumption  as  seeing  deteriorates,  and  Lee,  Holmes  and  Kerr  (1976)  have 
claimed  that  it  is  a  valid  assumption  at  least  in  the  absence  of 
turbulence.  Some  years  earlier,  Strohbehn  (1968)  commented  that  the 
Rayleigh  distribution  for  amplitude,  and  thus  the  negative  exponential 
distribution  for  intensity,  is  valid  for  "tropospheric  bevond-the- 
honzon  propagation  or  m  line-of -sight  propagation  when  the  turbulent 
medium  is  a  small  slab  and  the  receiver  is  far  from  the  slab."  This, 
however,  does  not  describe  the  situation  for  ground-based  astronomy. 
Certainly,  the  main  advantage  of  using  negative  exponential  statistics 
rather  than  lognormal  statistics  is  that  it  allows  one  to  more  easily 
achieve  analytic  solutions  to  many  of  the  problems  arising  in  astronomi¬ 
cal  imaging.  We  shall,  however  with  one  exception,  use  the  lognormal 
distribution  for  all  of  our  calculations  in  subsequent  chapters. 
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REVIEW  OF  SPECKLE  IMAGING 

At  this  poirt,  we  have  stated  that  the  atmosphere  degrades 
images  of  objects  viewed  through  it  by  distributing  a  point  source  such 
as  a  star  into  a  speckle  pattern,  and  we  have  determined  a  statistical 
model  for  the  intensity  of  an  optical  wave  which  has  passed  through  the 
atmosphere.  We  now  wish  to  discuss  the  previous  work  that  has  been 
done  on  the  problem  of  obtaining  true  images  of  objects  viewed  through 
the  atmospnere. 

This  problem  has  received  considerable  attention,  especially  in 
the  last  fifteen  years  since  the  advent  of  speckle  interferometry.  In 
the  remainder  of  this  chapter,  we  will  review  some  of  the  major 
advances  in  speckle  imaging.  In  section  3.1,  we  will  briefly  discuss 
short  exposure  images  followed  by  comments  on  conventional  long 
exposure  and  speckle  interferometry  in  sections  3.2  and  3.3.  Then 
sections  3.4  through  3.8  will  be  devoted  to  short  reviews  of  various 
techniques  used  to  extract  more  information  from  the  output  of  speckle 
interferometry  or  from  the  short  exposure  images  themselves.  There 
are,  of  course,  many  other  techniques  which  will  not  be  covered  in  this 
chapter,  and  we  refer  the  reader  to  Dainty  (1984)  or  to  Bates  (1982) 
for  a  more  thorough  review. 
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3.1  The  Short  Exposure  Image 
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A  short  exposure  or  "instantaneous"  image  refers  to  an  image 
for  which  the  fluctuations  of  the  atmosphere  may  be  considered  to  be 


jjfv; 

frozen;  that  is 

■ 

•  v*- 

of  a  second  .  VI 

*3 

m-z- 

systems  model  (’ 

i 

- 

i(x 

^  y 

where  i(x,y)  is 

object  intensit; 

i(x,y)  =  o( x, v )  **  t(x,y) 


(3.1) 


atmosphere/telescope  point  spread  function,  and  **  denotes  convolution. 

This  equation  may  be  equivalently  expressed  m  the  Fourier 
domain  bv 


*  - 


I(u, v)  a  0(u. v)T( U, v) 


(3.2) 


where  I(u,v)  and  Ci(u,v)  are  respectively  the  Fourier  transforms  of  the 
image  and  object  intensities  and  T(u,v)  is  the  combined 
atmosphere/telescope  transfer  function.  Since  we  are  dealing  with 
astronomical  imaging,  it  should  be  noted  here  that  although  Equations 
(3.1)  and  (3.2)  are  written  in  the  spatial  and  spatial  frequency 
domains,  it  is  more  appropriate  to  interpret  x  and  y  as  angles  of  arc 
and  u  and  v  as  angular  frequency,  that  is,  arc  sec 

Thus  our  problem  is,  given  i(x,y)  (or  equivalently,  I(u,v)),  to 
determine  o(x,y)  or  the  best  possible  image  of  o(x,y)  given  our  optical 
system.  In  the  short  exposure  case,  one  may  not  simply  use  the  point 
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Figure  3. 
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Linear  system  model  of  image  formation  through  the  atmos¬ 
phere.  —  (a)  Pictorial  description  of  image  formation; 
and  (h)  Block  diagram  of  the  linear  systems  model  of  image 
f  ormat i on . 
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spread  function  t(x.y)  or  the  transfer  function  T(u.v)  to  restore  the 
image  because  due  to  the  action  of  the  atmosphere,  t  (  x ,  y )  and  T  ( u ,  v ) 
are  random,  and  information  is  available  only  for  the  average 
quantities,  lt(x,y)>  or  <T(u,v)>,  This  would  imply  that  some  form  of 
averaging  of  the  short  exposure  images  is  necessary.  In  section  3.2, 
we  will  discuss  one  form  of  averaging,  the  long  exposure  image,  and  see 
that  it  also  has  serious  drawbacks  in  terms  of  loss  of  high  frequency 
information.  We  will  see  that  a  better  approach  is  the  shif t-and-add 
algorithm  of  section  3.8,  which  consists  of  averaging  properly 
registered  short  exposure  images  and  retains  information  out  to  the 
diffraction  limit. 


3.2  Long  Exposure  Imaging 

Conventional  long  exposure  imaging  may  be  regarded  as  a  sum  of 
a  series  of  short  exposures,  so  the  governing  equation  for  this  process 
in  the  Fourier  domain  is  as  follows 

<I( u , v) >  =  0(u,v)<T(u,v)>  (3.3) 

where  < . >  denotes  an  ensemble  average,  <I(u,v)>  is  the  Fourier 
transform  of  the  long  exposure  image,  and  <T(u,v)>  is  the  long  exposure 
transfer  function. 

Earlier  analyses  of  the  long  exposure  transfer  function  were 
performed  by  Hufnagel  and  Stanley  (1964)  and  Fried  (1966),  and  Fried 
has  shown  that  long  exposure  images  retain  substantially  less  high 
frequency  information  than  do  the  individual  short  exposures,  thus 
yielding  a  blurred  or  smoothed  result.  Intuitively,  we  may  see  this  bv 
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recalling  that  T(u.v)  is  random  due  to  atmospheric  effects  and  may  thus 
be  positive  or  negative  (or  perhaps  complex-valued)  at  high 
frequencies.  Summing  the  individual  T(u,v)’s,  as  in  Equation  (3.3). 
will  then  result  in  suppression  of  nigh  frequencies,  so  that  even  if 
<T(u,v)>  is  known,  one  cannot  adequately  reconstruct  o(x,y)  due  to  this 
lack  of  high  frequency  information  in  the  image  data. 

For  further  discussion  of  the  long  exposure  transfer  function 
and  determination  of  its  functional  form,  see  Fried  (1966),  or  Hufnagel 
and  Stanley  (1964)  or  the  review  articles  Dainty  (1984)  or  Roddier 
(1981). 

3.3  Speckle  Interferometry 

It  was  to  combat  this  loss  of  high  frequency  information  that 
Labeyrie  (19"0)  introduced  the  technique  which  became  known  as  speckle 
interferometry  (Gezarie,  Labeyrie,  and  Stachnik,  1972).  Where  long 
exposure  imaging  is  equivalent  to  summing  a  series  of  short  exposures 
and  thus  their  Fourier  transforms,  speckle  interferometry  consists  of 
the  addition  of  the  squared  magnitudes  of  the  short  exposure  image 
transforms  (i.e.,  spatial  power  spectra)  as  follows 

Gj(u,  v)  =  <  |I(u,v)  |2> 

-  |0(u,v) |2<iT(u,v)  12> 

■  $Q(u,v)< |T(u,v)  j2>  (3.4) 

where  $j(u,v)  and  4>q(u,v)  are  respectively  the  image  and  object  average 

2 

(spatial)  power  spectra  and  ].|  denotes  the  squared  magnitude  of  a 

i  2 

complex-valued  quantity.  This  averaging  of  |T(u,v)|  ,  which  is  always 
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a  non-negative  quantity,  prevents  the  loss  of  high  frequencies,  which 
was  the  main  drawback  of  long  exposure  imaging. 

1  ,  o 

Korff  (19’3)  has  derived  an  analytic  expression  for  <\7(v,v) 
based  on  lognormal  intensity  statistics,  which  agrees  closely  with 
Fried's  results  (Fried,  1966).  However,  this  expression  requires 
considerable  numerical  evaluation,  which  limits  its  usefulness.  Dainty 
(1984)  has  shown  that  a  good  approximation  for  <|T(u,v)i“>  is 

< | T(u , v) | 2>  »  | <T( u, v) >J  2  +  kTD(u,v) 

where  <T(u,v)>  is  the  long  exposure  transfer  function,  T_(u,v)  is  the 
diffraction-limited  optical  transfer  function  of  the  telescope,  and 
k  <  1  is  a  constant  inversely  related  to  the  number  of  speckles  in  the 
image. 

Now,  assuming  that  information  about  the  transfer  function 
<[T(u,v)j“>  is  available,  either  from  a  reference  star  (point  source) 
or  from  one  the  theoretical  derivations  mentioned  above,  <5>q(u,v)  may  be 
recovered  from  Equation  (3.4).  Inverse  transforming  yields  the  average 
spatial  autocorrelation 

<Rg(x,y)>  =  F_I{«>0(u,v)} 

*  <i(x,y)  68  i( x, y )> 

where  66  denotes  autocorrelation. 

The  limitation  of  this  method  is  that  one  obtains  only 
autocorrelation  or  power  spectrum  images  rather  than  true  images  of  the 
object.  Phase  information  is  lost,  and  as  is  noted  in  Oppenheim  and 
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Lim  (197u),  this  information  is  generally  more  important  than  amplitude 
information  for  signal  reconstruction.  Bates  (1982)  has  presented  a 
comprehensive  review  of  whac  has  become  known  as  the  "phase  problem," 
the  inability  (in  general)  to  uniquely  determine  an  object  o(x.y)  from 
its  power  spectrum.  In  certain  special  cases,  however,  such  as  a 
centre-symmetric  object,  the  Fourier  transform  is  purely  real  and 
determination  of  o(x,y)  is  possible.  Also,  other  useful  object 
information,  such  as  the  distance  between  two  point  sources  (e.g., 
double  stars)  or  an  estimate  of  the  spatial  extent  of  the  object,  nav 
be  recovered  from  the  autocorrelation  data. 

Since  the  introduction  of  speckle  interferometry  in  1970,  much 
effort  has  been  centered  on  methods  of  recovering  an  object  from  its 
autocorrelation  or  power  spectrum,  and  we  will  discuss  several  of  these 
algorithms  in  the  following  sections. 

3.4  Speckle  Holography 

The  technique  of  speckle  holography  (Bates,  Gough  and  Napier, 
1973;  Gough  and  Bates,  1974)  demonstrates  that  the  object  o(x,y)  may  be 
reconstructed  from  interferometry  data,  provided  an  unresolvable 
reference  object  (point  source)  such  as  a  star  is  present.  In  this 
case,  following  the  notation  of  Dainty  (1984),  the  object  may  be 
represented  as  a  sun  of  two  parts 

o(x,y)  »  6(x)<5(y)  +  o1(x-x1,v-y1) 

where  6(x)6(y)  is  the  Dirac  delta  function  denoting  the  reference 
object  and  o, (x-x^,y-y^)  centered  at  (x^.y^)  is  the  object  of  interest. 
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If  x,  N  3x  '2  and  v,  >  3v  /2.  where  x  and  v  are  the  extent  of  o.(x.v) 
I  o  '  1  '  o  o  '  o  1 

m  the  x  and  v  directions,  respectively,  then  the  spatial  autocorrela¬ 
tion  R  (x,y)  separates  into  three  distinct  parts 

R  (x,v)  =  o, (-( x+x. ) ,-( v+v. )  + 
o  i  1  ■  '  l 

[ 5( x )6( v )8S5( x ) 5( v )  +  o, (x. y)9@o. ( x. y )  ]  ♦ 

■i.  X 

o^((x-x^),(v-y^)). 

That  is,  the  central  component  of  R  (x,y)  consists  of  the  autocorre¬ 
lations  of  the  two  parts  of  the  object,  and  the  outer  components 
consist  of  their  cross-correlations.  One  of  these  outer  components 
will  be  the  correctly  oriented  object  of  interest  and  the  other  will  be 
a  180°  rotation,  so  the  object  can  be  reconstructed  within  this 
rotational  ambiguity. 

Practically,  the  object  may  not  always  meet  the  separation 
requirements  x,  >  3x  / 2  and  y,  >  3v  / 2  even  when  a  point  source  is 
present,  so  that  the  parts  of  RQ(x,y)  are  not  completely  distinct.  Liu 
and  Lohmann  (1973)  have  suggested  a  procedure  similar  to  that  of 
speckle  holography  in  which  they  utilize  the  fact  that  the  long 
exposure  image  does  contain  low  frequency  phase  information.  This 
information  is  incorporated  into  their  algorithm  by  using  the  long 
exposure  image  as  a  mask  to  select  the  correct  components  of  R  (x,y) 
and  eliminate  autocorrelation  terms. 
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Related  to  speckle  holography  is  the  speckle  masking  technique 


developed  by  Weigelt  (Wei gelt,  19?':  Weigel t  and  W'irnitzer,  1983),  who 


has  suggested  that  it  may  be  especially  useful  for  imaging  double 


stars.  Tne  major  difference  between  speckle  holography  and  speckle 


masking  is  that  in  the  former,  an  unresolvable  reference  object  is 


required  to  be  present  while  in  the  latter,  each  speckle  image  is 


preprocessed  nonLinearlv  to  create  a  synthetic  reference  object. 


Fo  implement  this  algorithm,  one  first  calculates  the  average 


image  triple  correlation 


<[ i(x, y) i( x-m^, y-m^) )00i( x, y ) > 


(3.5) 


where  (m  ,m  )  is  the  masking  vector,  e.g.,  in  the  case  of  imaging  a 
x  y 


double  star,  (m^.m  )  is  the  separation  which  may  be  obtained  from 


interferometry  data.  Weigelt  has  then  computed  a  correction  term 


( Weigel t  and  Wirnitzer,  1983),  which  allows  the  calculation  of  the 


object  triple  correlation 


<[o(x,y)o(x-m  ,  y-m  )]50o(x,y)> 
x  y 


from  Equation  (3.5).  The  object  o(x,y)  may  then  be  reconstructed 


provided  that  (m  ,m  )  was  properly  chosen  to  give 
x  y 


’o( x, y)o(x-mx, y-my) ]  *  5(x)5(y) 


3.6  Knox-Thomnson  Method 

Knox  and  Thompson  (1974)  have  proposed  a  method  that  involves 
the  autocorrelation  of  the  Fourier  transform  of  the  speckle  images 
rather  than  the  image  power  spectra.  The  algorithm  calculates 

<I(Ui,v1)I*fu,.v,)>  =  Ofu^,  )0*(u^,  v^XTfu^,  )T*(u0,  v0)> 

=  0(u,  v  )0*(  u+Au,  v+Av)<T(  u ,  v  )T*(  u+Au,  vi-lv  )  - 

where  Au  *  u0-u^  and  Av  *  v^-v^  are  small  compared  to  the  correlation 
length  of  I(u,v).  Denoting  the  phase  of  0(u,v)  by  4>(u,v)  and  the  phase 
of  T(u,v)  by  6(u,v),  we  may  write 

<I(u,,v, )I*(u,,v0)>  *  |0(u,v)  j jO(u+Au, v+Av)  j 

*  exp(  i  [  d>(  u ,  v )  -  6( u+Au . v+Av) ] ) 

*|T(u,v) ;  |T( u+Au, v+Av) I 

*exp( i[ 9(u, v)  -  9( u+Au, v+Av )]) . 

(3.6) 

Dainty  (1984)  has  shown  that  the  phase  difference  [ 9(u, v)-S( u+Au, v+Av ) ] 
=  0,  so  that  the  phase  of  the  right-hand  side  of  Equation  (3.6)  is 
approximately  [<Ku,v)  -  <t>( u+Au ,  v+Av)  ] .  Thus  we  may  generate  a  grid  of 
phase  differences  for  the  object  and  may  obtain  the  relative  phase  for 
each  point  in  the  Fourier  transform  of  the  object  by  summing  the  phase 
differences  between  it  and  an  arbitrary  reference,  say,  the  origin. 
Since  the  modulus  of  the  transform  is  available  from  the  interferometry 
data,  object  reconstruction  is  then  possible  (within  a  position 
ambiguity  since  only  relative  phase  may  be  computed). 


3. 7  Fienuo* s  Iterative  Algorithms 
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Fienup  (1978,  1979)  has  suggested  an  iterative  approach  to 


obtain  object  phase  information  from  the  modulus  of  the  Fourier 


transform  0(u,v)  ,  which  is  known  from  speckle  interferometry.  His 


approach,  called  the  error  reduction  method  (Figure  3.2)  because  the 


mean-squared  error  decreases  at  each  iteration,  is  a  modified  version 


of  the  Gerchberg-Saxton  algorithm  (Gerchberg  and  Saxton,  1972)  and 


consists  of  simply  changing  the  object  constraints  in  this  well-known 


method . 


Beginning  at  the  kLi  iteration  with  the  object  estimate 


°  A 

°k(x,y),  this  estimate  is  Fourier  transformed,  yielding  O^u.v)  = 
jO^(u, v) iexp[id^(u,v) ] .  The  Fourier  domain  constraint  consists  of 


replacing  |0  (u,v)'  with  the  known  modulus  |0(u,v)j.  The  quantity 


|0(u.v)  i exp [  if) ,(u,  v)  ]  is  then  inverse  Fourier  transformed  producing  the 


A  T 

image  o^(x,y),  which  is  then  forced  to  obey  the  spatial  domain 


constraints.  The  principal  spatial  domain  constraint  is 


non-negativity,  although  other  a  priori  information  about  the  object 


may  also  be  included,  e.g.,  we  know  the  object  autocorrelation  from 


speckle  interferometry  and  the  object  diameter  cannot  exceed  half  the 


extent  of  the  autocorrelation.  Thus  the  new  object  estimate 


Vi(*.y> 


o^(x,y),  constraints  satisfied 


constraints  not  satisfied 


Figure  3.2.  Fienup's  error  reduction  method. 


is  formed.  This  entire  procedure  may  be  start®:  “  a: 

v,(u,v)  or  a  random  i,(u,v)  if  a  better  estimate  is  una vai .a: . e . 

1  X 

In  practice,  Fienup  has  found  that  although  convergence  ;> 
initially  rapid,  it  soon  becomes  extremely  slow,  requiring  an 
impractical  number  of  iterations  for  a  good  reconstruction. 

To  combat  this  problem,  he  has  developed  the  input-output 
approach  (Figure  3.3)  in  which  the  new  object  estimate  is  a 
modification  of  the  previous  one 


(  x ,  y ) 


o,(x,y),  constraints  satisfied 


o,(x.v)  -  ao, ( x, y ) ,  constraints  not  satisfied 

K  K 


where  a  is  a  constant.  He  has  further  found  by  experimentation  that 
the  swiftest  convergence  is  achieved  by  periodically  varying  the  method 
of  forming  o.  . (x,v)  after  everv  few  iterations. 

°  Km  1  ' 


3.8  Shif t-and-Add 

Unlike  the  methods  previously  discussed,  the  shif t-and-add 
algorithm  and  the  related  method  which  will  be  examined  in  this  section 
do  not  attempt  to  reconstruct  the  object  from  its  power  spectrum  or 
autocorrelation  or  to  use  other  information  about  the  Fourier  modulus, 
.mage  separation,  or  image  extent  which  are  made  available  by  speckle 
interferometry.  Instead,  they  attempt  to  obtain  a  true  image  of  the 
object  directly  from  the  short  exposure  images. 

Shift-and-add  is  an  extension  of  a  method  previously  developed 
by  Lynds,  Worden  and  Harvey  (1976)  (or  Worden,  Lynds  and  Harvey,  1976), 
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wrv-  r»al izec  t r.at  eacn  s Deckle  of  a  short  exposure  image  is  itself  a 
distorted  irnao-1  of  toe  object.  Their  procedure  involves  the  identifi¬ 
cation  and  superposition  cf  the  brightest  speckles,  thus  creating  an 
improved  estimate  of  the  object. 

According  to  3ates  (19S2),  this  procedure  is  limited  by  the 
fact  that  many  celestial  objects  are  so  faint  that  a  typical  short 
exposure  image  contains  few  speckles,  and  he  and  Cady  (Bates  and  Cady, 
19SQ;  Cady  and  Bates,  1980)  have  extended  the  method  to  compensate  for 
this  problem.  In  their  approach,  each  short  exposure  is  shifted  so 
that  its  brightest  speckle  is  located  at  the  origin,  and  all  such 
shifted  images  are  summed,  thus  the  name  shif t-and-add .  Assuming  that 
this  processing  is  also  carried  out  for  a  reference  star,  the 
shi f t-and-add  image  of  the  object  of  interest  may  be  deconvolved  using 
the  shif t-and-add  image  of  the  reference  star  as  a  point  spread 
funcr Lon  . 

The  beauty  of  this  method  is  in  its  simplicity  and  ease  of 
implementation  as  compared  to  the  methods  previously  discussed,  and 
according  to  Bates  (1982),  it  may  be  digitally  implemented  in  real 
tine.  Further  improvement  in  the  output  image  has  been  obtained 
Oates,  1982)  by  a  method  called  adjusted  shift-and-add ,  in  which  each 
pixel  is  multiplied  by  the  value  of  the  brightest  pixel  before  the 
summing  operation,  and  Bates  and  Robinson  (1982)  have  devised  a 
slightly  more  complicated  version  of  shift-and-add,  which  they  have 


found  useful  m  ultrasonic  imaging. 


We,  however,  will  restrict  ourselves  to  me  original ,  simple 


version  of  the  algorithm  and  will  present  a  theoretical  analysis  of  the 
method  and  various  simulation  results  in  the  following  chapters. 


CHAPTER  4 

ANALYSIS  OF  THE  SHI FT- AND- ADD  ALGORITHM 
FOR  LOGNORMAL  INTENSITY  STATISTICS 

The  simplicity  of  the  shi f t-and-add  algorithm  described  in 
section  3.S  leads  us  to  ask  why  this  method  forms  an  improved  estimate 
of  the  object.  This  question  has  previously  been  addressed  in  some 
detail  by  Hunt,  Fright  and  Bates  (1983)  for  negative  exponential 
statistics.  Here  we  shall  present  an  analysis  (which  draws  heavily  on 
the  earlier  one  for  notation  and  modelling)  using  the  lognormal 
intensity  statistics  which  were  justified  in  Chapter  2. 

To  facilitate  our  discussion,  we  will  now  adopt  the  standard 
mathematical  model  of  the  speckle  process.  We  are  changing  notation 
from  that  of  section  3.1  to  be  more  consistent  with  existing 
shi f t-and-add  literature.  This  model  assumes  that  image  exposure  time 
is  so  short  that  the  atmospheric  amplitude  and  phase  variations  are 
essentially  frozen  and  also  makes  the  assumption  of  isoplanaticity , 
i.e.,  the  assumption  that  the  atmospheric  point  spread  function  is 
shift  invariant.  With  no  loss  of  generality  and  to  simplify  notation, 
we  shall  restrict  our  analysis  to  one  dimension.  We  denote  the  mth 
speckle  image  of  the  object  f(x)  by 

s  (x)  *  h _(x)*f(x)  +  c  ( x )  (4.1) 

mm  m 
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where  m  is  a  tine  index,  h  (x)  is  the  short  exposure  speckle  ooint 

m 

SDread  function,  *  indicates  convolution,  and  c  (x)  is  a  general 
‘  m 

contamination  term  that  includes  all  other  degradations,  i.e., 

recording  noise,  system  nonlinearities,  photon  noise. 

Within  this  framework,  the  shift-and-add  process  is  easily 

described.  It  consists  of  finding  the  maximum  value  of  each  spatial 

image  and  the  spatial  coordinate  at  which  it  occurs,  translating  the 

image  bv  £  so  that  the  maximum  value  is  now  located  at  the  origin  of 
m  ° 

coordinates,  and  summing  all  such  translated  images.  For  M  speckle 
images,  the  shift-and-add  result  s(x)  is  expressed  by 

M 

s(x)*l/M)s(x+£)  (4.2) 

m  m' 
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With  this  model  established,  we  proceed  to  explore  our  original 

question:  why  is  s(x)  a  better  estimate  of  f(x)  than  are  any  of  the 

individual  frames  s  (x)? 

m 

We  first  consider  this  question  qualitatively.  According  to 
Bates  and  Cady  (1980),  the  brightest  portion  of  each  speckle  image 
s  (x)  is  likely  to  be  a  distorted  version  of  the  brightest  portion  of 
f(x).  Thus  because  of  the  linear  nature  of  the  degradation  model  and 
of  the  shift-and-add  process,  superposition  of  distorted  versions  of 
the  brightest  portions  of  f(x)  necessarily  implies  the  correctly 
registered  superposition  of  all  other  parts  of  f(x). 


•  « 
f 


4.1  Derivation  of  the  Point  Spread  Function 
Now  we  perform  a  more  rigorous  analysis  of  the  shift-and-add 
method  by  considering  the  case  of  a  single  point  source  (Dirac  delta 
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function)  object  and  finding  an  overall  point  spread  function  for  the 
degradation  plus  shif t-and-add  processing.  We  will  at  this  point 
ignore  the  contamination  term  c  (x)  so  that  Equation  (4.1)  becomes 


s  (x) 
m 


=  h  (x) . 
m 


Thus  the  overall  point  spread  function  h(x)  is 


M  M 

h(x)  =  s ( x )  =  1/M  £  s  (x  +  K  )  =  1/M  l  h  (x  +  E  ) 

m  m  m  m 


m=l 


m=l 


as  is  obvious  from  the  linearity  of  the  method. 

Following  the  derivation  in  Hunt,  Fright  and  Bates  (1983),  we 
introduce  a  change  of  notation 


c 

m 


(x) 


S  (x  + 
m 


r 

'’m 
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i.e.,  the  am(x)  are  simply  the  speckle  images  translated  so  that  the 
maximum  value  is  found  at  the  origin.  This  leads  to  the  shi f t-and-add 
result 

M 

s( x )  =  1/M  l  a  (X). 
m=l  m 

At  each  fixed  x,  this  is  simply  a  sum  of  M  random  variables,  where  we 

note  that  for  | x|  very  small,  correlation  in  the  images  necessitates 

that  o  (x)  -  0  (0),  while  for  |x|  larger  than  the  correlation  length  of 
mm 

the  speckle  process,  am(x)  and  0(0)  are  effectively  independent. 
Therefore,  for  large  M, 

M 

lim  s(x)  *  lim  1/M  7  a  (x)  ->  E[o  (x)|c  (0)1 

LLwrt  ^  IT1  IT1 
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where  £[•]  denotes  expectation,  and  we  will  assume  M  large  enough  that 


£ 


lim  s(x)  =  E[o  (x) ic  (0)] 


M-k» 


m  m 


(4.3: 


in  the  following  calculations.  We  also  define  for  future  use  the 
quantity  s 


max 


M 


s  =  1/M  I  o  (0) 
max  m 

m=l 


or  for  M  large,  we  mav  assume  s  =  Elc  (0)1. 

°  1  max  m 


Again,  following  the  notation  of  Hunt,  Fright  and  Bates  (1983), 


we  let 


h  -  V01 


i2  -  vx)- 


As  we  are  assuming  lognormal  intensity  statistics,  1^  is  governed  by 
the  probability  density  function 


Pdi) 
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exp 


-[iogi1-u1r 
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and  similarly  for  I^.  In  this  expression,  log  denotes  the  natural 

2 

logarithm  (base  e),  and  u  ^  and  are  respectively  the  mean  and 
variance  of  the  normally  distributed  random  variable  log  I ^ . 

To  calculate  the  expectation  in  Equation  (4.3),  we  need  the 
joint  density  function  for  I.  and  I,,  which  we  will  assume  to  be 
jointly  lognormal: 


p(IrI2)  = 


r0>l\'ui\  /iogu-u,  /  u 

,1  °i  /  l  °2  / ' r 
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where  o  is  the  correlation  coefficient  for  the  associated  normal  random 
variables  log  1^  and  log  I7.  We  have  made  such  an  assumption  because 
it  is  necessary  in  order  to  achieve  analytic  results  although  we  are 
well  aware  that  it  is  not  generally  valid  to  do  so  simply  given  that 
the  variables  1^  and  I ^  are  both  marginally  lognormal.  Also  we  note 
that  the  x  dependence  of  p  and  of  is  not  explicitly  stated  in 
Equation  (4.4),  and  we  will  comment  further  on  this  later. 

Then  with  our  previous  assumptions,  we  have 

,  PCIj-1!’ 

p(I2'Il)  *  ~pa7) 
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and  we  may  write  the  desired  expectation 


E[ljl,  ]  «  I,p(lJl,)dI. 
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This  integral  is  evaluated  in  Appendix  B,  giving 
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To  give  some  intuitive  feel  for  the  expression,  we  will 
evaluate  it  for  some  extreme  values  of  the  parameters.  Also  needec  for 
this  evaluation  are  the  following  relations  (see  Appendix  A)  for  anv 
lognormal  random  variable  I: 


<I> 


E[I]  - 


D4C2/2 

e 


and 


■> 

=  'I>2(eC  -1) 


•-here  as  before 


V  =  E[ log  I  ] 

2 

=  var[ log  I ) . 

First  we  consider  the  case  where  and  1^  are  completely  uncorreiated . 

Recalling  that  I,  =  a  (0)  and  I_  *  c  (x),  this  is  simplv  the  case  of  !x' 
i  m  L  m 

>>  0.  This  corresponds  to  the  parameter  value  c  *  0.  Thus  Equation 
(4.6)  becomes 

u^c;/2 

E[  1 2 1 1 ^ 1  *  e  ‘  ‘  -  <I2> 

which  is  exactly  what  one  intuitively  expects. 


Next  we  consider  the  case  of  I0  and  1^  perfectly  correlated, 
i.e.,  X;  -  0.  This  corresponds  to  the  parameter  value  p  =  1.  For 
this  case,  it  is  also  completely  reasonable  to  make  the  further 
assumption  that  U1  2  and  0“  =  clj.  Then  evaluation  of  Equation  b\ 
v  leids 


1J  '  "1 


which  accords  exactly  with  intuition. 

At  this  point  we  observe  that  Equation  (4.6)  is  entirely 
specified  by  parameters  associated  with  the  normal  distributions  of  log 
I,  and  log  I?.  We  shall  assume  that  the  speckle  images  are  spatially 
stationary,  or  at  least  wide-sense  stationary,  an  assumption  actually 
valid  only  near  tne  origin  (or  only  locally),  which  implies  a  constant 
mean  and  variance 


<!,>  *  <I^,>  *  <I> 
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which  in  turn  implies 


'1  *  yn  *  u 


Having  made  this  assumption,  our  further  analysis  will  be  strictly 
valid  oniv  in  the  region  of  x  -  0.  Since  we  are  particularlv 


*  '.‘‘'.V.V.V,*.*  V.V  / 
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interested  in  this  region--to  see  how  closely  s(x)  resembles  tne 
original  object  -(x) — our  results  should  yield  the  desired  information. 


SH 


Returning  to  our  notation  of  d  (0)  and  "  (x)  rather  tnan  ..  and 


we  now  nave 


E'~  =.  x)  i -  (0)]  =  c  (0)  exp[t,-ou+c"(  1-cT) -'f 
'  m  '  m  m 


=  a  (0)<I>1_0  expt^-^—o2' 

m 


(4.7) 


and  our  next  step  is  to  relate  C  =  0(x),  the  correlation  coefficient  of 

log  b  (0)  and  log  c  (x),  to  r(x),  the  correlation  coefficient  of  a  (0) 

°  m  °  m  n 

anc  ~  (x).  Bv  definition, 

m 

cov [ I _ , I , 3 

r  (  x  )  =  — - — - - 

-1^1, 

A.  A. 

V1: 

so  we  will  now  determine  the  quantity  E[  I ^  I ] .  We  make  the  same 
assumption  that  I,  and  I0  are  jointly  lognormally  distributed  as  was 
made  in  the  calculation  of  Thus 
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This  integral  is  evaluated  in  Appendix  B.  giving 
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Making  our  same  stationaritv  assumptions,  which  imply  that 
we  obtain 

c(x)c2  , 

,  ^  e  -i 

r(  x)  =  - ? - 
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c  *  e(x)  *  l/c2  log [  1  +  r(x)(ec^  -  1)] 

Now 

i  0,  Jl 

cf  =  <!>“[ e-^  -  1 ] , 
so 

c2  =  log[l  +  c'/CLi4-] 

and  after  some  algebraic  manipulation 

0-  o(x)  ,  i°g£W.--.^g<l^ 
logR(o]  -  log<I>‘ 

where  R(x)  denotes  autocorrelation.  Assuming  ergodicity,  R(x)  =  R^ 
the  spatial  autocorrelation  of  the  speckle  image. 

Substitution  in  Equation  (4.7)  and  further  algebraic 


simplification  yields 
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and  since  we  are  implying  that  M  ->  00  by  using  the  expectation,  we  have 
the  desired  result 


lir.  s(x) 
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The  obvious  question  at  this  point  is  how  to  characterize  R  (x) 

in  order  to  give  this  complicated  expression  some  meaning.  There 

appear  to  be  three  options  available.  First,  one  may  actually  measure 

R  (x)  from  the  speckle  images.  We  had  no  real  data  available  but  have 

included  plots  of  s  (x)  and  the  corresponding  R  (x)  taken  from  images 
m  s 

simulated  by  the  method  of  section  5.1  in  Figures  4.1  and  4.2. 

n 

Korff  (1973)  has  developed  an  expression  for  <|T(f)|“>  which  is 

the  Fourier  transform  of  <Rs(x)>.  Calculating  <|T(f),  >  and  inverse 

Fourier  transforming  is  our  second  option.  Unfortunately,  evaluation 
,  n 

of  <T(f)|->  requires  extensive  numerical  computation  and  specific 


parameters  for  the  atmosphere  and  optical  system,  so  we  shall  consider 
our  third  possibility. 

This  third  option  for  characterizing  S  (x)  is  an  approximation 
to  Korff's  result  developed  by  Dainty  (1975).  Dainty  has  claimed  that 
Korff's  result  is  in  "broad  agreement"  with  that  obtained  by  the 
simpler  model 

< : T ( f ) i 2 >  »  ;<T(f)>;2  +  kTD(f) 

i  2 

where  <T(f)>,~  is  the  squared  modulus  of  the  long  exposure  transfer 
function,  T^(f)  is  the  diffraction-limited  transfer  function,  and  k  <  1 
is  a  constant  depending  on  atmospheric  and  telescope  parameters.  As  in 
Hunt,  Fright  and  Bates  (1983),  we  inverse  Fourier  transform  this 
equation,  obtaining 

<Rs(x)>  =  1 C  x )  +  l<Ka(x) 

where  i(x)  is  the  long  exposure  component  of  Rs(x),  i.e.,  the 
autocorrelation  of  the  seeing  disc;  a(x)  is  the  diffraction-limited 
component  of  R  (x),  i.e.,  the  Airy  disc  of  the  telescope;  and  K  is  a 
normalizing  factor  required  to  force  <Rs(0)>  =  <!““>. 
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Then  from  Equation  (4.8), 


1  lit  s(x)  =  s 
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log<I“>  -  log<I>“ 


(4.9) 


The  first  term  of  this  expression  is  a  constant,  so  diffraction-limited 
behavior  will  occur  in  the  second  term  depending  on  the  relative 
magnitudes  of  l(x)  and  kKa(x);  that  is,  if  the  diffraction-limited 
kKa(x)  is  large  enough  compared  to  the  broad  smooth  function  l(x). 

This  discussion,  together  with  Equation  (4.9),  completes  our 
characterization  of  the  point  spread  function  of  the  shif t-and-add 
process . 

4.2  Analysis  for  Two  Point  Sources 
The  next  degree  of  complication  in  this  procedure  is  obviously 
to  consider  an  object  consisting  of  two  point  sources.  The  behavior  of 
the  shif t-and-add  process  is  unchanged  although  now  the  possibility  of 
ghost  peaks  exists,  and  it  is  the  probability  of  such  an  occurrence 
that  will  be  the  major  focus  of  this  section. 

Vie  first  define  our  object 

f(x)  *  a^5(x-x^)  +  a25(x-X2>  (4.10) 

where  a^  >  a2*  Then  clearly  each  speckle  image  takes  the  form 


(*.. in 


s  (x)  =  ah  (x-x  )  *  a.h  (x-xn), 
m  I  m  i  I  m  - 

where  we  are  again  ignoring  the  contamination  term.  Shi ft-and-add 
processing  will  allow  objects  to  be  resolved,  each  witn  a  profile  of 
the  form  of  Equation  (4.8),  if  the  distance  ! x0-x, .  between  the  two 
point  objects  is  greater  than  the  correlation  length  of  the 
diffraction-limited  component  of  the  speckle  image.  At  this  point  we 
define  some  notation 


w  =  a ,,  /  a  1 

<I,>  =  mean  value  of  a, h  (x-x,) 
1  1ml 

<10>  =  mean  value  of  a.h  (x-x~) 
-  -mu 


=  variance  of  a.h  (x-x.) 

1  m  1 

=  variance  of  a_h  (x-x,.) 

2  m  2 

2  2 


where  of  course  we  are  concentrating  our  attention  on  the  central 
portion  of  each  speckle  pattern  ^^(x-x^)  where  the  speckle  image  may 
be  treated  as  a  stationary  random  process  with  constant  mean.  With 
this  notation  established,  we^'hnay  use  Equation  (4.8)  to  define  a 
composite  shi ft-and-add  profile 


For  the  shift-and-add  process  to  be  carried  out  correctly,  the  maximum 

value  of  s  (x)  should  be  the  maximum  of  a, h  (x-x.).  There  are  two  wavs 
m  1  m  1 

that  this  may  be  prevented  from  occurring.  First,  the  sum  of  any  two 

speckles  from  a. h  (x-x,)  and  a-,h  (x-x-)  may  exceed  the  sum  of  the 
1  m  1  2  m  2 

maximum  of  a.h  (x-x.)  and  any  speckle  from  a-h  (x-x-).  Intuitively, 
i  m  i  2  m  2 

one  would  not  expect  this  to  occur  in  any  systematic  fashion,  since 

this  would  be  completely  unrelated  to  the  maximum  of  a-h  (x-x-).  Thus 

2  m  2 

if  this  error  occurred  repeatedly,  it  would  result  in  a  general 

randomness  in  s(x)  rather  than  in  a  false  peak.  The  second  way  that  an 

incorrect  maximum  may  be  chosen  will  result  in  a  ghost  peak  if  it 

occurs  repeatedly.  This  happens  when  the  maximum  of  a-h  (x-x-)  plus 

2  m  2 

anv  speckle  from  a.h  (x-x.)  exceeds  the  sum  of  the  maximum  of 
i  m  1 

a.,  h  (x-x,)  and  any  speckle  from  a-h  (x-x_).  The  probability  with  which 
l  ill  i  L  m  z 

this  occurs  will  determine  the  relative  magnitude  between  correct  and 
false  peaks. 

Therefore  we  will  calculate  this  probability  using  the  notation 
of  Hunt,  Fright  and  Bates  (1983).  We  let 

z  =  maximum  of  a.h  (x-x.) 

I  m  1 

<I2>  a2 


wz  =  maximum  of  a-h  (x-x-) 

2  m  2 

r.v.(l)  =  a  random  variable  with  probability  density 
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P(I7)  - 


I ,  C  ,  V' 4TT 
1  1 


exp  1 


-[loglj-uJ 

20f 


=  a  speckle  from  a.h  (x-x,), 
lml 


and  r.v.(2),  likewise. 


V.'e  Chen  define 

u  =  z+r . v . ( 2 ) 
v  s  wz+r . v . ( 1) 


The  probability  that  the  correct  maximum  of  s^(x)  will  be  chosen  is 

P  [u>v]  and  the  probabilitv  that  the  maximum  of  a.h  (x-x.)  will  be 

/  m  L 

chosen  instead  is  P  [v>u]  =  1  -  P  [u>v].  (We  are  ignoring  here  the 

possibility  of  two  speckles  adding  together  to  exceed  either  u  or  v.) 

Now,  for  each  speckle  image,  both  the  magnitudes  and  positions 
of  z  and  wz  are  fixed.  Also,  given  our  assumption  that  the  two  point- 
source  objects  are  separated  by  more  than  the  correlation  length  of  the 
diffraction-limited  component  of  the  speckle  images,  we  may  consider 
r.v.(l)  and  r.v.(2)  to  be  independent  but  conditioned  on  z.  Thus 

p(u. v|z)  =  p(u jz)p(v \z) 


which  from  simple  probability  is 


p(u, v jz) 


2*0,02  (U"-Hv_w:) 


exp 


~  [  log  (u-;)  -U-,] ' 


2c! 


Vie  will  now  simplify  this  expression  by  recalling 

<I.>  *  mean  value  of  a.h  (x-x.) 

1  lml 

<I->  »  mean  value  of  a.h  (x-x.)  »«<'!'• 

2  2  m  2  1 


] log (v-wz) -u  ^  ]  *" 


CT "  =  variance  of  a,h  (>:-x  ) 

1 .  »  m  i 

1  O  O 

CT"  =  variance  of  a“h  (x-x~)  =  w"<I, 
X  ^  m  a.  « 


and  for  lognormal  random  variables 
<lx>  = 


-i  ?  w  1 

=  <  I .  >  -  [  e  -  1],  likewise  for  ■' 

ll  i  ' 


This  leads  to 


~>  -> 
C2  =  C1 


U-,  =  U,  +  log  W, 


so  we  set 


I 


and 


°1  =  C2  =  C 

U1  =  U 

=  U  +  log  w 


i 


J-[log(”)-u]'-[log(v-w-:)-u]l 


p(u’v'2)  =rb-(uT7y(V-zT  exp 


2c 


Then  the  probability  integral  of  interest  is 


cov 

„r  .  ,  i  ff  i  r-riog(v-wt)-u]2\  i  j 

P(V>U>  -  ^2  J  —  -P  { - 

wzz 


- 1 - >dudv . 

7n~  J 


We  first  consider  the  case  of  w  =  1,  for  which  the  integral  is  analytic¬ 


ally  evaluable.  Then 


and  we  make  me  change  of  variable 

s  =  log' v-z) 
t  =  log(u-z) 


to  obtain 


P  [v>u] 


-CD 


rt-ur 
”  -> 

e  dtds 


and  using  the  evaluation  of  the  normal  probability  integral  from 
Abranovitz  and  Stegun  (1970) 


P  [v>u]  = 


2a 


CO  - 

1  [  2a2 

n:  e  l 1 

v  ?7T  j 


er  f(— x)  ]ds 

<Jy  2 


Vie  now  let 


P  * 


s-u 


and  the  integral  becomes 

00 

1  t 

P  [  v>u  ]  »  -  I  dpe 


2  r  2 

_p  +  — ~  I  dpe  P  erfp. 


2»'Tr  J  J 

»QO  —CD 

The  second  integral  is  equal  to  zero  because  its  integrand  is  odd,  so 
our  result  is  P  [v>u]  *  1/2,  which  agrees  entirely  with  our  intuition. 


The  obvious  question  to  address  at  this  point  is  how  fast  this 
probability  decreases  as  w  decreases  from  1.  Unfortunately,  for  values 


of  w  o  trier  than  1,  the  integral  is  ana  :  v t : c a  1 1  ••  intractable  so  we 
resort  to  numerical  integration.  Numerical  integration  requires 
calculation  of  this  probability  in  trie  form  P  ;  v  u ;  =  1  -  P  ;u  u; 
a vo ic  1  o  car 1 1  b.Ti s  of  necai :  ve  numbers.  Thus  we  are  10  calculate 


- exp 

J  J  v-w: 


- i iOg 1 v-wc 


[loglHlcct-J' 


d  v  au 


and  making  the  usual  change  of  variable 
s  =  log(v-wz) 


,u*c. 

t  =  iog(— -). 


the  integral  becomes 


p  !  u>v  ]  «  _L_  '  [  ex 

2~o“  J  J 


logfwe1* (l-w)c) 

I  exp{h?L}  ”p  1 


Evaluation  of  the  s  integral  yields 


?  [  u>v  ]  *  - —  j  exp 

2a»/2Tr  ‘ 


mi 


1  *  erf  (  Ut 

av3 


2o»'5r  l  2o“  >  (  oO  J 


We  let  p  *  — -  and  obtain 
o/2 
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INTENTIONAL  LEFT  OUT 


To  avoid  evaluation  or  3  double  integral,  we  have  used  the 
following  expression  from  Abranowitz  and  Stegur.  (  1970}  as  an 
approximation  for  erf  x: 


where 


=  1  -  (a 


+  a  ..t  -“-at  +  a . 


le  X  ef x 


t  = 


1 


l*px 

I  e(  x )  j  _<  1.5  x  10 
p  =  .3275911 
a  =  .254829592 
a0  =  -.284496736 


a3  =  1.421413741 
a4  =  -1.453152027 
a5  =  1.061405429 


The  integrand  of  Equation  (4.12)  is  plotted  in  Figure  4.3.  and  the 
integration  results  for  three  sets  of  parameters  are  summarized  in 
Tables  4.1  through  4.3. 

The  probabilities  of  selecting  the  wrong  maximum  seem 
unrealistically  low,  although  they  may  indeed  be  realistic  given  the 
relatively  large  separation  between  the  original  point  sources  that 
assumed  in  the  modelling.  However,  we  do  in  part  attribute  the  low 
probabilities  to  the  absence  of  noise  and  nonlinearities  m  the 
derivation  and  possibly  to  our  assumption  of  constant  means  and 
variances  <I,>,  < 1 n > ,  c.  ,  c,  .  This  essentially  is  an  assumption 


Table  4.1.  Probability  of  selecting  wrong  naximum: 

average  parameter  values.  —  I^>  =  II. 6; 

=  27.1;  and  z  =  243.55.  <!,>  and 

for  this  table  are  the  average  values  of 
the  parameters  (100  sets)  supplied  by  Jon 
Freeman . 


P  (  v  >  u ) 


% 
r  V 

>V 

W 

Second  Max 

wz 

Simpson's  Rule 

7-Point 

Gaussian 

Quadrature 

r 

i 

243.55 

.500 

.500 

.  99 

241.11 

.356 

.356 

I; 

.95 

231.37 

.052 

.052 

t 

.9 

219.19 

.0035 

.0035 

1 

.85 

207.02 

.0003 

.0003 

X 

.8 

194. S4 

.00003 

. 00004 

.75 

182.66 

. 000004 

.00001 

.5 


121.78 


.000001 


.00001 


1 


Table  4.2.  Probability  of  selecting  the  wrong  maximum: 

parameter  values  from  a  saturated  frame.  — 
<I.>  =  15.72;  o j ,  =  37.55;  and  z  =  255. 


w 

Second  Max 

wz 

P  (v  >  u) 

7-Point 

Gaussian 

Simpson's  Rule  Quadrature 

1 

255.00 

.500 

.500 

.  99 

252.45 

.375 

.375 

.95 

242.25 

.075 

.  072 

.9 

299.59 

.0063 

.0063 

.35 

216.75 

.0006 

.00059 

.8 

204 . 00 

. 000064 

. 0000~2 

.75 

191.25 

. 0000089 

.000017 

.5 

127. 50 

.000001 

. 0000099 

i  a  b i e  4.3. 


Probability  of  selecting  the  wrong  maximum: 
parameter  values  from  an  unsaturated  frame.  — 
< I ^ >  =  9.55;  =  16.44;  and  z  =  218.00. 


w 

Second  Max 

v/z 

P  (v  >  u) 

7-Point 

Gaussian 

Simpson's  Rule  Quadrature 

1 

218.00 

.500 

.500 

.  99 

215.82 

.334 

.334 

.95 

207.17 

.0342 

.0342 

.9 

196.20 

.00166 

.00166 

.85 

185.30 

.0001 

.0001 

.8 

174.40 

.00001 

. 00002 

.  75 

163.50 

. 000002 

.00001 

.  5 

109.00 

.000001 

.00001 
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st at ionar i t v  which  is  in  actuality  valid  only  in  a  small  region  about 
the  maxima  z  and  wz. 

In  spite  of  these  unrea  lis'.ical  lv  low  probabilities,  we  may- 
still  draw  the  important  conclusion  that  the  correct  maximum  will  be 
chosen  in  the  majority  -.f  f-ames  and  the  shif t-and-add  process  carried 
out  correctly  even  when  the  original  point  sources  do  not  differ 
greatly  in  magnitude. 


4.3  Analysis  of  the  Extended  Object  Case 

We  now  address  the  case  of  most  practical  interest,  the  case  of 

an  extended  object  rather  than  a  set  of  isolated  point  sources.  We 

then  rewrite  Equation  (4.1)  for  this  case 

00 

S  (x)  =  [  f(x-x  )h  (x. )dx  +  c  ( x )  (4.13) 

m  j  i  m  I  l  m 

*00 

as  a  convolution  integral,  again  following  the  notation  of  Hunt,  Fright 

and  Bates  (1983).  Ignoring  the  contamination  term,  s  (x)  is,  for  fixed 

m 

x,  an  integral  of  lognormal  random  variables  h  (x. )  weighted  bv  the 

mi 

object  itself,  f(x-x^).  Digitally,  of  course,  this  becones  a  weighted 
sum  of  lognormal  random  variables,  or  just  a  sum  of  nonidentical ly 
distributed  lognormal  random  variables. 

Thus  the  relevant  question  is:  what  are  the  statistics  of  sums 
of  lognormal  random  variables?  For  the  case  of  N  independent  random 
variables,  the  sum  tends  to  become  Gaussian  distributed  as  N  approaches 
infinity  bv  the  Central  Limit  Theorem.  Barakat  (1976)  has  investigated 
this  case,  with  the  added  assumption  that  the  random  variables  are 
identically  distributed.  His  approach,  as  summarized  below,  is  a 


Also  due  co  independence,  the  characteristic  function  of  Z  is 


Bn 


tz  (t)  =  j  eXZkZ  P(zk)dzk 

^  -oc 

-  (it)n<(Xk-E[Xk]) 

=  i  n 

n=0  (varX)  n! 


Fhus,  with  representing  the  i  central  moment  of  Z, 

.3.  4 


rz  AL  3  4  _5/7  N 

4  (t)  =  d-%-  -J-37^  +  - ”  +  0(N  )}‘  ’ 

2  -N  60JN  2404N’“ 


so  taking  the-  logarithm  off  both  sides  and  using  the  expansion 


log(l-s)  =  l  (-l)k+1  4-.  -1  <  s  <  1 

k=l 


a  s  -  1/2  S2  +  1/3  S3, 


we  obtain 


2  iV-t  V3  4  3/2 , 

log  4  ( t )  =  -  4  - T~ 1  +  0(N  } 

2  1  6c“N  ' “  240  N 


. .2/2  if  ft  _  ^ 

4  (t)  =  e-C  '“(1  -  —^75  +  0(N  *)) 
2  6Ni/2 


v,  =  - —  is  the  coefficient  of  skewness. 

1  ,,3/2 


Fourier  transforming  yields 


d.,  (  2 )  =  — —  e  2  ^  ( 1  +  — ■  l~  h,(z)  +  0(N  )) 


where  h^(z)  denotes  the  Hermite  polynomial  z  -3z.  Thus  for  large 

finite  N.  sums  of  independent  lognormal  random  variables  converge  to  a 

-1/2 

Gaussian  probability  density  function  only  as  N  .  This  result 
should  also  hold  qualitatively  for  non-identical ly  distributed  random 
variables  as  long  as  the  Central  Limit  Theorem  conditions  are  met  (see 
Sect  ion  2.1.1). 

Our  problem  is  complicated  further  by  the  fact  that  our 
lognormal  random  variables,  h  (x),  arc  weakly  correlated  as  was 
illustrated  in  Figure  (4.2).  We  see  from  this  figure  that  the 
diffraction-limited  portion  of  the  autocorrelation  is  quite  narrow; 
indeed  the  width  of  the  base  of  the  narrow  peak  is  only  twice  the 
extent  of  a  speckle.  Since  Barakat ' s  result  is  based  on  use  of  the 
Central  Limit  Theorem,  we  assume  that  it  would  also  apply  (perhaps  with 
even  slower  convergence)  for  random  variables  which  are  only  weakly 
correlated.  This  is  due  to  the  fact  that  variations  of  the  Central 
Limit  Theorem  do  exist  for  dependent  random  variables  (Loeve,  1950; 
SerfLing,  1968;  Hoeffding  and  Robbins,  1948).  According  to  Serfling 
(1968),  the  assumptions  required  for  these  Central  Limit  Theorems  are 
in  practice,  not  amenable  to  proof.  It  is  further  stated  that  many 
experimenters  feel  that  the  Gaussian  approximation  is  valid  for 
stationary  processes  which  are  observed  for  time  (in  our  case  spatial) 
intervals,  which  are  long  in  comparison  to  the  correlation  length  of 
the  process. 

We  have  assumed  that  our  images  are  stationary,  at  least 
locally,  and  have  noted  that  the  diffraction-limited  correlation  length 
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Implicit  in  this  calculation  is  the  assumption  of  stationanty  of  the 
speckle  point  spread  function  and  the  integrated  speckle  images,  which 
we  continue  to  emphasize  is  valid  only  near  the  origin. 

Inverse  Fourier  transforming  Equation  (4.16)  yields  the 
autocorrelation  function  R  (x)  in  terms  of  the  previously  discussed 
(section  4.2)  R  (x)  =  P*h(x)  and  the  autocorrelation  of  the  object 

Re(x)  -  Rf(x)*Rh(x).  (4.17) 

As  we  noted,  R, (x)  may  be  decomposed  into  a  low  frequency  and  a 
diffraction-limited  component 

R,  (x)  2  l(x)  +  kKa(x) 
h 

and  we  now  suppose  that  the  object  power  spectrum  may  also  be 
decomposed  into  low  frequency  and  diffraction-limited  components 

I F( u) | 2  =  | Fl(u) | 2  +  i Fd(u) 1 2 

or  in  the  soatial  domain 

' 

Rf(x)  «  R1(x)  +  Rd(x). 

Substituting  in  Equation  (4.17)  gives  us 

R  (x)  *  l(x)*[R,(x)+R.(x)]  +  kKa(  x)*R,  (  x)  +  kKa(  x)*R  .( x) . 
e  Id  x  d 

(4.18) 

Thus  if  the  object  contains  a  spatial  frequency  component  out  to  the 
diffraction  limit,  the  profile  of  s(x)  in  Equation  (4.15)  may  also 
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* 

i 


I 
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contain  a  dif f raction- limited  component  provided  that  the  magnitude  of 

kKa(x)*R^(x)  is  large  enough  relative  to  the  remainder  of  Equation 

(-.18).  Addition  of  the  correlation  component  due  to  the  contamination 

term  c  (x)  mav  further  obscure  the  diffraction-limited  information, 
m 

Now  we  return  to  the  case  of  a  large  object  f(x)  and  assume  a 

Gaussian  distribution  for  the  speckles  of  s  (x).  The  analysis  for  this 

m 

case  has  been  presented  by  Hunt,  Fright  and  Bates  (1983).  As  in 
section  4.1,  the  conditional  expectation  is  calculated; 

however,  the  conditional  density  function  now  takes  the  form 


P(I2|IL)  = 


aJ2^(l-T)2 


exp 


l  2o  (i-r"( 


(I2-<Ie»-r(l 


r<!e>,]:) 


(4.19) 


where  <Ig>  and  o^  are  respectively  the  mean  and  variance  of  the 
integrated  speckle  pattern  and 


r  =  r(x) 


Re(x) 

- *r~ 

C  ~ 


e 

is  the  normalized  correlation  of  the  integrated  speckle  pattern.  We 
further  note  that  r(x)  is  a  zero  mean  correlation  coefficient  with 
r(x)  <_  1,  and  that  the  usual  stationarity  properties  are  assumed. 
Evaluation  of  E [  I ^  i  I -^  ]  gives 


lim 

M-h» 


s  (x)  = 
e 


<1  > 

e 


e(x) 


[s  -  <1  >] 
max  e 


(4.20) 


e 

where  Re(x)  is  characterized  as  in  Equation  (4.18)  but  with  the 

2 

quantity  <Ig>  subtracted  because  the  correlation  is  required  to  be 
zero  mean .  Thus 
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,(x)  =  l(x)*{R1(x)+Rd(x)]+kKa(x)*R1(x). 


■kKa(x)*R  ,(x)-<I 
d  e 


(4.21) 


and  if  desired,  one  r,ay  also  include  the  autocorrelation  of  the 

contamination  term  R^fx).  As  before,  we  will  obtain  diffraction- 

limited  resolution  if  the  magnitude  of  kKa(x)*Rd(  x)  is  large  enough 

relative  to  the  other  components  of  Equation  (4.21).  In  both  Equations 

(4.1  =  )  and  (4.21),  1( x)*[ R ^( x)  +  Rrf( x )  ]  and  kKa(x)*R^(x)  are  smooth 

functions  and  the  useful  information  kKa(x)*R  ,(x)  will  thus  be  found 

d 

riding  on  a  broad  background. 

So,  as  noted  by  Hunt,  Fright  and  Bates  (1983),  if  Rg(x) 
possesses  a  large  enough  diffraction-limited  component  and  f(x) 
contains  a  strong  point  source,  the  linearity  of  the  shif t-and-add 
process  will  reconstruct  the  extended  object  (on  a  background).  Also, 
a  very  strong  point  source  should  ensure  that  the  maximum  value  is 
correctly  chosen  for  each  sm(x).  Such  sources  may  be  expected  to  occur 
often  in  natural  imagery,  e.g.,  an  unresolvable  star  near  an 
astronomical  object  of  interest  or  a  glint  on  the  wing  of  an  airplane. 

In  Chapter  6,  we  present  supporting  simulation  results  for 
shifr-and-add  processing  of  a  simulated  astronomical  object  with  a 
nearby  "star"  and  demonstrate  that  further  improvement  may  be  obtained 
by  Wiener  filtering  of  the  shif t-and-add  result. 

4.4  Miscellaneous  Results 

In  the  remainder  of  this  chapter  we  will  examine  other  results 
on  the  shift-and-add  algorithm.  In  section  4.4.1,  we  will  analyze  the 
question  of  convergence  which  is  of  considerable  interest.  Following 


this,  we  will  discuss  some  of  the  work  of  Bagnuolo  (1984)  in  which  he 
examines  the  effect  of  nonisoplanat ic it y  (i.e.,  lack  of  shift 
invariance  over  the  entire  image)  on  the  shi f t-and-add  method. 

4.4.1  Convergence  of  the  Shi ft-and-Add  Algorithm 

In  sections  4.1  through  4.3  we  have  analyzed  the  shi f t-and-add 
method  and  found  that  it  may  produce  an  image  containing 
diffraction-limited  information.  In  all  of  this  analysis,  it  was 
effectively  assumed  that  an  infinite  number  of  images  was  summed  which 
is,  of  course,  not  feasible  in  practice.  Thus  we  are  led  to  inquire 
just  how  fast  the  shif t-and-add  algorithm  converges — does  it  require 
five  or  five  thousand  summations'* 

To  answer  this  question,  we  recall  that  for  each  fixed  x 
shift-and-add  consists  of  summing  M  lognormally  distributed  random 
variables,  where  M  is  the  number  of  frames  processed.  These  random 
variables  are  independent,  since  they  occur  in  different  frames;  and 
since  we  are  considering  a  fixed  spatial  position,  we  will  also  assume 
that  they  are  identically  distributed.  Hence,  Barakat's  (1976)  results 
(as  discussed  extensively  in  section  4.3)  apply  directly.  That  is,  the 

—  l/o 

sum  of  random  variables  approaches  its  limiting  distribution  as  M 

We  stress,  as  does  Barakat,  that  this  conclusion  is  valid  only  for 

large  M.  We  nave  observed  in  our  simulations  that  convergence  is  nore 
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rapid  for  small  M,  say  <  10.  However,  this  would  be  true  even  if  M 
convergence  is  valid. 

We  may  see  this  same  result  another  way,  by  considering  a 
"noise-to-signal  ratio"  for  the  sum  of  independent  identically 
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distributed  lognormal  random  variables  at  fixed  spatial  position  x. 
That  is,  we  are  constructing  a  ratio  of  the  spread  about  the  mean  to 
the  mean  of  the  sum.  For  convergence,  we  wish  this  ratio  to  decrease 
to  zero.  Following  the  notation  of  section  4.1  and  ignoring 
contamination  effects,  we  have  the  sum 


M 


s(x)  *  I  <Mx) 

m=l  n 

which  has  mean  and  variance 


E[s(x)]  =  ME[o  (x)]  =  M<o  (x)> 
m  m 


2,  C 


var  [  s(  x)  ]  =  Mvar[cr  (x)  ]  =  M  <o  (x)>“[e  -lj 


where  =  var [log  c  (x)].  Then  our  noise-to-signal  ratio  for  the  s 
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The  numerator  of  this  quantity  is  constant,  so  the  noise-to-signal 
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ratio  decreases  as  M  ,  which  agrees  with  our  prior  conclusion. 

Thus  we  have  determined  that  convergence  of  the  shif t-and-add 
algorithm  is  slow;  however,  in  our  simulations  (Chapter  6),  notable 
improvement  has  been  achieved  with  sequences  of  only  twenty  frames.  An 
important  conclusion  suggested  by  this  slow  convergence  is  that  it  may 
well  be  advantageous  to  average  fewer  frames  and  achieve  further 
deblurring  by  other  techniques;  for  example,  Wiener  filtering,  rather 
than  simply  continuing  to  increase  M. 
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4.4.2  Effect  of  Nonisoplanat icit y  on  Shi ft-anc-Add 

For  the  entire  preceding  analysis,  we  have  made  the  assumption 
of  isop lana t lei t y .  That  is,  we  have  assumed  tnat  the  image  was  formed 
by  passage  of  wavefronts  from  the  entire  object  through  one  statistic¬ 
ally  homogeneous  portion  of  the  atmosphere,  giving  us  a  shift-invariant 
point  spread  function.  According  to  Dainty  (1*375),  it  has  been  found 
experimentally  that  this  assumption  is  invalid  for  star  fields  or 
larger  objects. 

In  this  case,  Equation  (4.13)  becomes 


s  (x)  = 
m 


f( x ' )h  (x; x  '  )dx ’ 
m 


(4.22) 


where  we  have  ignored  the  contamination  term  c  (x).  The  intuitive 

n 


approach  is  to  decompose  hm(x;x’)  into  a  spatial  sum  of  shift  invariant 


point  spread  functions 


h  (x)  •  I  h  (x-x  ) 
m  mj  j 


where  K  is  the  number  of  isoplanatic  patches  through  which  the  object 
wavefronts  pass.  One  would  then  make  convergence  arguments  allowing 
the  change  of  order  of  summation  and  integration  in  (4.22),  yielding 


s  (x)  "Is  .(x-x  . ) . 

m  i-m-i  j  ' 


j= 1 


Then  shi f t-and-add  gives 


1  M 

S(x)  "  M  I, sm(x) 
m=l 


1  M  K 

rr  I  J  s  .(  x-x  . ) . 

Mm=l  j=l  mJ  J 
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At  this  point  the  analysis  becomes  Intractable  because  K  is  a  random 
variable  depending  on  the  step  m.  Therefore,  another  approach  should 
be  taken  in  order  to  achieve  any  results. 

Bagnuolo  (1984)  has  addressed  this  problem  for  the  specific 
case  of  the  shif t-and-algorithm  applied  to  imaging  of  double  stars.  In 
an  earlier  paper  (Bagnuolo,  1982),  he  has  presented  a  slightly  modified 
shif t-and-add  algorithm  and  a  method  for  calculating  the  intensities  of 
the  background,  the  two  valid  peaks  representing  the  double  star,  and 
me  accompanying  ghost  peak.  Ratios  of  these  intensities  are  then  used 
to  compute  the  true  intensity  ratio  for  the  two  components  of  the 
couble  star.  In  his  later  publication  (Bagnuolo,  1984),  he  has  shown 
mat  his  method  is  invariant  to  the  degree  of  isoplanaticit y .  This  is 
by  no  means  a  general  theoretical  result,  but  it  nay  indicate  that 
nonisoplanaticity  is  not  a  serious  problem  in  shift-and-add  processing. 
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CHAPTER  5 

COMPUTER  SIMULATION  OF  TURBULENCE  DEGRADATIONS 

In  order  to  test  the  effectiveness  of  the  shift-and-add 
algorithm,  especially  for  the  case  of  extended  object  imaging,  we  next 
performed  simulations  of  the  image  degradations  produced  by  atmospheric 
turbulence.  Series  of  point  spread  functions  have  been  computed,  both 
for  an  algorith.  developed  by  McGlamerv  (1976)  that  considers  only 
phase  perturbation;:  of  the  optical  wave  and  for  a  modified  version  of 
this  algorithm,  which  also  includes  amplitude  perturbations.  These 
point  spread  functions  were  then  convolved  with  a  simulated  object,  and 
the  resulting  series  of  degraded  object  frames  were  subjected  to 
shift-and-add  processing  and  other  restoration  techniques,  as  discussed 
in  Chapter  6.  In  this  chapter,  we  will  examine  the  two  algorithms  for 
computing  point  spread  functions  that  we  mentioned  above. 

5.1  Point  Spread  Function  Simulation 
for  Phase  Perturbations  Alone 

The  first  algorithm  that  we  used  for  point  spread  function 
simulation  was  one  developed  by  McGlamerv  (1976).  In  this  algorithm, 
he  has  considered  only  phase  perturbations  of  the  optical  wavefront 
since  he  claims  that  these  are  the  dominant  factor  in  image  degrada¬ 
tion.  He  also  assumes,  as  is  common  (Strohbehn,  1968)  that  the 
Kolmogorov  spectrum 
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is  a  valid  description  of  the  spatial  power  spectrum  of  the  phase  at 
the  receiver,  where  C  is  a  constant  related  to  turbulence  strength  and 
wavelength  and  f  denotes  spatial  frequency.  With  these  preliminary 
assumptions  noted,  we  now  summarize  the  algorithm. 

Step  one  consists  of  generating  a  complex  array  of  Gaussian 
random  numbers  that  are  to  represent  the  spatial  frequency  domain  for 
the  phase  of  the  wavefront.  (Phase  is  assumed  to  be  Gaussian 
distributed,  as  discussed  in  Chapter  2,  which  necessarily  implies  that 
it  is  also  Gaussian  in  the  frequency  domain.)  In  step  two,  we  multiply 
the  array  by  f-1^^,  the  square  root  of  the  Kolmogorov  spectrum.  This 
factor  introduces  correlation  into  the  point  spread  function,  for 
without  it,  one  produces  a  random  noise  field  rather  than  a  speckle 
pattern.  The  third  step  consists  of  Fourier  transforming  the  array  and 
separating  the  result  into  its  real  and  imaginary  parts.  Each  of  these 
arrays  represents  a  map  of  phase  at  the  pupil  of  the  optical  system, 
which  we  denote  by  <f>( u , v )  .  In  step  four,  we  form  the  path  length 
difference  map 


d( u , v) 


<b  (u,v)A 


Tn 


where  1  denotes  wavelength,  and  we  note  that  this  step  is  unnecessary 
if  one  is  considering  monochromatic  illumination.  The  intensity  point 
spread  function  is  calculated  in  step  five  as  follows 

h(x,y,A)  *  |F  {p  ( u ,  v)exp(  i2ird(  u  ,  v )  /  a)  *  J2 


denotes  the 


where  p(u.v)  is  the  pupil  function  of  the  telescope,  ;  •  | 
squared  naanitude  of  a  complex  quantity,  and  Ft’;  refers  to  the  Fourier 
transfer”.  The  final  step  six  is  not  computed  in  the  monochromatic 
case.  It  involves  forming  the  polychromatic  point  spread  function 

h(x,y)  =  k  ^  F(X)h(x,y,X) 

A 

where  k  is  a  normalizing  factor  required  so  that  the  value  of  the 
transfer  function  at  zero  frequency  is  unity,  and  F(X)  is  a  weighting 
function  which  depends  on  spectral  distributions  of  source  power, 
atmospheric  transmittance,  and  the  image  sensor.  McGlamery  (1976)  also 
discusses  other  factors  necessary  for  scaling  in  computations  of  the 
polychromatic  point  spread  function. 

In  our  simulations,  we  omitted  steps  four  and  six,  and  we 
varied  the  standard  deviation  of  the  Gaussian  iterates  in  step  one  in 
order  to  vary  the  strength  of  the  turbulence,  rather  than  changing  the 
constant  C  of  the  Kolmogorov  spectrum.  Furthermore,  we  ignored  the 
effects  of  the  telescope  pupil  and  set  p(u,v)  identically  equal  to 
unity.  Examples  of  point  spread  functions  simulated  in  this  manner  are 
shown  in  Figure  5.1.  It  should  be  observed  that  these  point  spread 
functions  are  of  greater  spatial  extent  than  they  appear  in  this 
figure.  The  outer  edges  do  not  appear  because  of  a  lack  of  dynamic 
range  in  the  photographic  process. 


For  the  sake  of  completeness  and  also  because  amplitude 

perturbations  are  important  in  the  assumption  of  lognormal  intensity 

statistics  (see  Chapter  2),  we  modified  McGlamerv's  algorithm 

(McGlamery,  1976)  to  include  the  effect  of  amplitude  perturbations. 

According  to  Strohbehn  (1968),  the  log-amplitude  spatial  power  spectrum 
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is  also  proportional  to  the  Kolmogorov  spectrum  f  '  ,  and  we  know 

from  Chapter  2  that  log-amplitude  is  also  governed  by  a  Gaussian 
probability  distribution.  Therefore,  simulation  of  amplitude  effects 
proceeds  in  much  the  same  manner  as  the  simulation  of  phase  effects  in 
section  5.1. 

Step  one  consists  of  generating  a  complex  array  of  Gaussian 
ran  numbers  that  represents  the  spatial  frequency  domain  for 
log-amplitude.  In  step  two,  the  array  is  multiplied  by  f  ^  the 
square  root  of  the  Kolmogorov  spectrum,  in  order  to  introduce  correla¬ 
tion.  We  then  Fourier  transform  the  array  and  separate  the  result  into 
its  real  and  imaginary  parts  for  step  three.  Each  of  these  arrays 
represents  a  map  of  log-amplitude,  denoted  logA,  at  the  telescope 
pupil.  Using  the  phase  map  <t>(  u,v),  which  is  generated  by  the  method  of 
section  5.1  and  assuming  monochromatic  light,  we  complete  the  fourth 
and  final  step  of  the  algorithm,  the  formation  of  the  intensity  point 
spread  function 


o 


i  ,  2 

h  (  x ,  y )  =  ,F -’p(u,  v)exp(  logA  +  i$(u,v)}  “ 

2 

where  as  before,  ?{•}  denotes  Fourier  transform,  j • indicates  squared 
magnitude  of  a  complex-valued  quantity,  and  p(u,v)  is  the  pupil 
function  of  the  telescope  which  was  set  identically  equal  to  unity  in 
all  of  our  simulations. 

Examples  of  the  point  spread  functions  simulated  by  this 
procedure  are  presented  in  Figure  5.2.  The  phase  variations  in  this 
figure  are  statistically  the  same  as  those  of  Figure  5.1,  so  any 
difference  is  due  to  amplitude  effects.  As  we  see  from  Figure  5.2,  the 
mam  effect  that  the  amplitude  perturbations  have  is  a  general  blurring 
of  the  point  spread  function,  and  therefore,  of  any  image  convolved 
with  it.  This  will  be  more  clearly  seen  in  Figure  6.3  of  the  next 
chapter . 


5.3  Effect  of  Setting  One  Pixel 
Equal  to  One  Speckle 

The  simulation  of  point  spread  functions  discussed  in  sections 
5.2  and  5.3  results  in  speckle  patterns  in  which  each  speckle  is 
represented  by  a  single  pixel.  In  reality  each  speckle  has  a  spatial 
extent  on  the  order  of  the  size  of  the  Airy  disc  of  the  telescope 
(Dainty,  1975),  so  what  one  is  actually  simulating  by  the  above 
algorithms  is  integrated  or  averaged  speckle.  As  we  have  assumed 
throughout  that  the  intensity  at  each  point  of  the  speckle  pattern  is 
distributed  lognormallv,  we  have  a  sum  of  partially  correlated 
lognormal  random  variables.  This  is  exactly  the  situation  of  section 


d.A,  where  the  weighting  function  is  equal  to  one.  As  in  this  previous 
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Figure  5.2.  Point  spread  functions  including  both  phase  and  amplitude 
degradations.  —  Phase  is  statistically  identical  to  that 
of  Figures  5.1a  and  5.1b. 

a.  Phase  corresponds  to  that  of  Figure  5.1a  with  severe  amplitude 
degradation . 

b.  Phase  corresponds  to  that  of  Figure  5.1b  with  less  severe  ampli¬ 
tude  degradation. 
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Figure  5.2a.  Phase  corresponds  to  that  of  Figure  5.1a  with  severe 
amplitude  degradation. 


rigure  5.2b.  Phase  corresponds  to  that  of  Figure  5.1b  with  less  severe 
amplitude  degradation. 
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analysis,  we  expect  that  this  integrated  or  averaged  speckle  will  be 
approximately  Gaussian  distributed  when  the  number  of  random  variables 
averaged  is  large  (Barakat,  1976),  but  will  be  better  described  by  a 
lognormal  distribution  when  fewer  random  variables  are  averaged 
(Mitchell,  1968).  As  we  are  averaging  over  an  area  the  size  of  the 
Airy  disc  in  this  case,  we  would  expect  relatively  few  pixels  to  be 
averaged,  implying  that  a  lognormal  distribution  is  appropriate.  In 
our  simulations,  telescope  effects  were  ignored,  which  is  equivalent  to 
assuming  an  infinite  pupil  and  thus  a  point  source  Airy  disc.  For  this 
case,  one  pixel  actually  represents  one  speckle;  but  in  more  realistic 
cases,  statistics  of  integrated  speckle  rather  than  simple  speckle 
should  be  used  in  simulations. 


CHAPTER  6 


RESTORATION  OF  SIMULATED  IMAGES 

In  chis  chapter  we  shall  attempt  to  restore  images  degraded  by 
atmospheric  turbulence  using  shift-and-add  processing  alone  and  then 
combining  shift-and-add  processing  and  Wiener  filtering  (actually 
pseudo-inverse  filtering  in  the  absence  of  noise)  in  order  to  confirm 
our  previously  obtained  theoretical  results.  Restoration  will  be 
possible  only  if  the  derivations  of  Chapter  A  are  correct,  that  is, 
only  if  diffraction-limited  information  is  indeed  preserved  by 
shift-and-add.  Otherwise,  residual  blur  will  be  apparent  in  the  final 
results  because  of  loss  of  high  frequency  information. 

Using  the  methods  described  in  Chapter  5,  we  have  generated 
series  of  twenty  point  spread  functions  (psf’s)  for  varying  degrees  of 
turbulence  and  have  convolved  these  psf's  with  an  undegraded  image.  We 
emphasize  that  in  these  simulations,  no  other  degradations  (e.g., 
photon  noise,  sensor  noise)  have  been  introduced,  although  they  are 
certainly  present  in  actual  imagery. 

In  a  previous  study  (Hunt,  Morgan  and  West,  1983),  we  have 
described  the  simulation  and  restoration  of  degraded  images  using  a 
simulated  silhouette  of  an  airplane  as  the  original  image.  However, 
Bates  pointed  out  that  the  use  of  such  a  large  object  with  a  small  psf 
is  a  contradiction  of  reality.  In  order  to  simulate  a  more  realistic 
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situation,  we  have  used  for  our  undegraded  image  a  small  circle  at  gray 
level  110,  witr.  a  dark  (gray  level  0)  and  a  bright  (gray  level  192) 
spot  m  the  interior.  The  circle  is  of  approximately  the  same  spatial 
extent  as  the  smallest  psf  considered,  and  we  added  the  interior  spots 
in  order  to  test  the  ability  of  our  processing  to  restore  detail  other 
than  the  object  outline.  The  undegraded  image  is  intended  to  represent 
a  fairly  large  astronomical  object,  and  therefore,  we  have  also  inclu¬ 
ded  a  "star"  (point  source)  at  gray  level  255  in  the  frame  in  order  to 
allow  us  to  estimate  the  overall  psf  of  atmospheric  degradation  plus 
shif t-and-add  processing.  This  psf  is  used  in  further  deblurring  after 
shif t-and-add .  Our  undegraded  image  is  pictured  in  Figure  6.1;  a 
series  of  degraded  images  with  only  phase  perturbations  (corresponding 
to  the  psf's  of  Figure  5.1)  is  shown  m  Figure  6.2;  and  both  phase  and 
amplitude  variations  (corresponding  to  the  psf's  of  Figure  5.2)  are 
included  in  the  degraded  images  of  Figure  6.3. 

6.1  Shif t-and-Add  Processing 

With  series  of  degraded  images  at  our  disposal,  our  next  step 
is  restoration,  which  will  first  be  attempted  by  shif t-and-add 
processing  alone.  Figure  6.1  depicts  the  results  of  shif t-and-add  of 
twenty  images  (corresponding  to  Figure  6.2),  which  have  been  distorted 
by  phase  variations  alone,  and  the  shif t-and-add  results  for  both  phase 
and  amplitude  degradations  (corresponding  to  Figure  6.3)  are  shown  in 
Figure  6.5.  We  note  that  there  is  comparable  improvement  in  the  case 
of  both  phase  and  amplitude  corruptions  relative  to  the  phase 


Figure  6.2.  Obiect  blurred  by  convolution  with  psi  s  containing  pnas 
degradations  only. 

a.  Blur  produced  by  psf  of  Figure  5.1a. 

b.  Blur  produced  by  psf  of  Figure  5.1b. 


igure  6.3.  Object  degraded  by  convolution  with  psf's  containing  botn 
phase  and  amplitude  degradations. 

Blur  produced  by  psf  of  Figure  5.2a. 

Blur  produced  by  psf  of  Figure  5.2b. 
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degradation  only  case;  and  just  as  one  would  expect,  the  restoration  is 
better  for  milder  turbulence. 

It  is  worthwhile  at  this  point  to  note  that  although  the  gray 
level  of  the  "star"  or  point  source  is  higher  than  that  of  the  bright 
interior  spot  in  the  object,  the  bright  spot  is  larger  in  area,  and  the 
maximum  value  chosen  from  each  frame  by  shift-and-add  is  a  speckle 
generated  by  the  bright  spot.  This  observation  implies  that 
shift-and-add  processing  may  well  lead  to  a  substantially  improved 
estimate  of  the  object  as  long  as  there  is  a  relatively  small  bright 
area,  not  necessarily  a  point  source,  present  in  the  object.  Thus 
shift-and-add  may  be  useful  processing  for  a  larger  class  of  natural 
imagery  than  simply  images  of  astronomical  objects.  One  would  expect 
the  registration  of  the  images  to  be  less  accurate  in  such  cases, 
resulting  in  a  more  blurred  shift-and-add  image;  however,  as  seen  in 
the  following  sections,  the  overall  point  spread  function  of  the  blur 
may  be  estimated  or  measured  and  then  removed  by  Wiener  filtering. 

We  nave  further  noted  confirmation  of  our  estimates  of  the 

convergence  of  this  algorithm  as  discussed  in  section  4.4.  In  a  prior 

report  (Hunt,  Morgan  and  West,  1983),  we  noted  a  comparison  of 

shift-and-add  results  for  twenty  and  fifty  frames  in  which  we  observed 

no  visible  improvement  in  the  fifty-frame  result  over  the  twenty-frame 

result.  This  accords  well  with  our  determination  that  the  algorithm 
-1/2 

converges  as  M  (where  we  recall  that  M  is  the  number  of  frames 
averaged ) . 
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Finally,  we  note  chat  while  there  is  considerable  improvement 
in  the  shi f t-and-add  results  over  a  single  short  exposure  frame  or  over 
a  simulated  long  exposure  image  (Figure  6.6),  there  is  still  a  sizable 
amount  of  blur  present.  The  following  sections  are  devoted  to  removal 
of  this  residual  blur. 

6.2  Determination  of  the  Overall  Point  Spread 
Function  for  Atmospheric  Degradations 
Plus  Snif t-and-Add  Processing 

In  order  to  apply  further  restoration  techniques  to  the 
shif t-and-add  image,  we  must  measure  or  estimate  the  point  spread 
function  (psf)  for  the  entire  process.  If  a  point  source  such  as  our 
simulated  "star"  is  included  in  the  original  object,  the  psf  may  simply 
be  extracted  from  its  shift-and-add  image.  This  is  the  approach  we 
have  used  in  our  restorations;  however,  in  cases  of  severe  degradation, 
subtraction  of  background  blur  from  the  psf  is  necessary. 

We  also  used  a  Gaussian  least-squares  fit  to  the  point  spread 
function  in  some  restorations.  This  is  because,  in  the  limit  of  large 
M,  we  expect  the  point  spread  function  to  approach  a  Dirac  delta 
function  (or  very  narrow  function)  riding  on  a  broad,  smooth  function 
as  in  Equations  (4.18)  or  (4.21)  and  the  following  discussion.  Thus  we 
are  assuming  a  Gaussian  form  for  this  smooth  function.  This  may  be  a 
useful  approach  for  estimation  of  the  psf  when  an  auxiliary  point 
source  is  not  present  in  the  object  if  one  can  find  a  valid  method  of 
determining  the  variance  of  the  Gaussian. 

If  there  is  no  point  source  present  in  the  object,  it  may  also 
be  possible  to  estimate  the  psf  from  the  edge  response  of  the  object 
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Comparison  of  shi f r-and-add  and  simulated  long 
Images.  —  Both  were  generated  from  a  series  of 
degraded  bv  phase  perturbations  alone. 

a.  Shift-and-add. 

b.  Long  exposure. 


9 


'igure  6.6. 


* 

* 

£ 


I 


exposu 
20  imag 


S8 


(Tatian,  1965;  Tescher  and  Andrews,  1972).  This  type  of  procedure  will 
work  best  if  the  object  is  fairly  large,  containing  large  sharp  edges, 
a  fact  which  will  limit  its  usefulness  in  astronomy.  However,  in  this 
case  a  nearby  scar  is  often  present,  as  was  assumed  in  our  simulations. 


6.3  Wiener  Filter  Restoration 

We  now  assume  that  we  have  obtained  our  desired  point  spread 
function,  whether  measured  or  estimated,  and  proceed  with  the 
deblurring.  In  our  previous  work  (Hunt,  Morgan  and  West,  1983),  we 
used  both  the  Wiener  and  Cannon  filters  (Andrews  and  Hunt,  1977)  but 
subjectively  determined  that  better  results  were  achieved  with  the 
Wiener  filter.  Therefore,  we  shall  use  the  Wiener  filter  in  our 
present  restorations  and  shall  employ  the  following  form  of  the  filter 


W(u, v)  = 


_ H* (u,v) _ 

iH(u,v)i2  +  On”/$£(u,v) 


where  H(u,v)  is  the  Fourier  transform  of  the  psf,  *  denotes  complex 

conjugation,  |.|  indicates  the  modulus  of  a  complex-valued  function, 

.  2 

Tj(u.v)  is  the  object  power  spectrum,  and  o  is  noise  variance.  We 

*■> 

have  not  actually  measured  or  estimated  a  value  for  c"  but  have  simplv 

n  ■ 

considered  a  constant  which  is  varied  to  achieve  the  best  restoration. 
Also,  we  used  the  actual  <t’^(u,v)  which  we  determined  from  the 
undegraded  image  although  in  practice,  this  is  not  known  and  must  be 
estimated  from  the  degraded  image.  Often,  one  would  have  some 
information  such  as  the  expected  object  shape  that  will  help  in  these 


estimations . 
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Now,  letting  S(u,v)  represent  the  Fourier  transform  of  the 
shi i t-and-add  image,  our  estimate  f(x.v)  of  the  object  is 

f(x,v)  =  F  ^  i S( u, v)W( u, v)} 

where  as  usual,  F  ^{ • }  denotes  the  inverse  Fourier  transform.  In  order 
to  achieve  good  results  in  cases  of  severe  degradation,  i.e.,  minimize 
"ringing,"  we  were  forced  to  smooth  the  sharp  edges  between  the  images 
and  the  zero  backgrounds  which  were  added  as  padding  for  the  Fourier 
transforms.  We  found  this  to  be  necessary  both  in  the  shif t-and-add 
image  and  the  measured  psf ,  and  it  was  accomplished  by  multiplying  the 
images  by  wide  Gaussians. 

These  restorations  have  been  performed  for  the  shif t-and-add 
images  depicted  in  Figures  6.4  and  6.5,  and  the  results  are  shown  in 
Figures  6.7  and  6.8.  From  these  figures  we  can  see  that  the  Wiener 
flit  r  has  done  an  excellent  job  of  removing  blur  left  in  the 
snift-and-add  images — the  outlines  of  the  object  are  sharp,  and  the 
interior  spots  have  been  well-reconstructed.  It  is  important  to  note 
here  that  this  ability  to  reconstruct  the  original  object  may  be 
regarded  as  confirmation  of  the  theoretical  results  of  Chapter  4 — 
diffraction-limited  information  is  indeed  preserved  by  shif t-and-add, 
at  least  in  the  absence  of  noise. 

We  have  also  performed  Wiener  filtering  in  a  few  cases  with  a 
Gaussian  least-squares  fit  to  the  measured  point  spread  function  (as 
discussed  m  section  6.2),  and  one  such  outcome  is  shown  in  Figure  6.9. 
Some  blur  still  remains  in  these  images,  and  we  attribute  this  to  the 


Figure  6.7.  Wiener  filter  restorations  of  the  shi f t-and-add  images 
depicted  in  Figure  6.4.  —  Original  image  was  degraded  by 
phase  perturb at  ions  alone. 
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igure  6.8. 


Wiener  filter  restorations  of  the  shi f t-and-add  images  of 
Figure  6.5.  --  Original  image  was  degraded  by  both  phase 
and  amplitude  perturbations. 

Wiener  filter  restoration  of  the  image  in  Figure  6.5a. 

Wiener  filter  restoration  of  the  image  in  Figure  6.5b. 
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fact  that  only  twenty  images  were  processed  by  shif t-and-add .  As  noted 
m  section  6.2,  this  Gaussian  approximation  should  improve  as  the 
number  of  frames  processed  increases.  Even  with  this  remaining  blur, 
the  interior  spots  are  clearly  visible,  and  we  have  observed  some 
improvement  over  shi ft-and-aad  alone. 

In  practice,  one  would  not  expect  such  good  results  as  obtained 
in  Figures  6.7  and  6.8  because  of  noise,  lack,  of  information  about  the 
exact  form  of  the  object  power  spectrum  ^(x,y),  possible  "ghosts"  in 
the  shi f t-and-add  image  and  possible  problems  in  obtaining  the  overall 
psf;  nevertheless.  Figures  6.7  through  6.9  indicate  that  shif t-and-add 
processing  plus  Wiener  filtering  shows  promise  for  restoring  images 
degraded  by  atmospheric  turbulence  or  other  random  media. 


v* 


CHAPTER  7 


CONCLUSION 


In  this  dissertation,  the  task  was  undertaken  to  analyze  the 
shift-and-add  method  in  order  to  provide  a  better  understanding  of  why 
it  creates  improved  versions  of  images  degraded  by  atmospheric 
turbulence. 

We  first  provided  the  necessary  statistical  framework  and  image 
model  for  this  analysis  along  with  a  brief  summary  of  other  imaging 
methods  employed  in  the  presence  of  turbulence.  Then,  in  Chapter  4,  we 
demonstrated  that  the  result  of  shift-and-add  processing  of  a  series  of 
short  exposure  images  may  contain  diffraction-limited  information,  both 
in  the  point  source  object  and  in  the  extended  object  cases.  U:e  also 
showed  that  the  probability  of  error  in  such  processing  is  relatively 
low,  although  the  results  were  obtained  for  a  noise-free  case. 

Further,  we  presented  a  convergence  rate  for  this  algorithm  which  is 
initially  fairly  rapid,  but  which  slows  considerably  as  the  number  of 
images  processed  increases.  This  result  suggested  the  restoration 
method  applied  in  Chapter  6:  use  shift-and-add  to  average  relatively 
few  images  and  achieve  further  improvement  by  another  method. 

In  support  of  this  analysis,  we  presented  simulation  results 
which  demonstrated  the  effectiveness  of  the  shift-and-add  method, 
especially  when  combined  with  Wiener  filtering.  When  one  considers 
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that  these  results  were  obtained  by  processing  only  twenty  inages,  the 
performance  of  the  algorithm  becomes  even  more  impressive. 

As  a  further  note  in  support  of  shi f t-ana-processing ,  our 
previous  work  (Hunt,  Morgan  and  West,  1983)  indicated  that  it  is  of 
value  in  reducing  the  error  m  features  extracted  from  turbulent 
images.  The  specific  feature  extraction  algorithm  considered  was 
invariant  moments  (Hu,  1962),  and  in  almost  all  cases  we  found  that  the 
moments  of  the  shift-and-add  processed  images  were  lower  in  error  than 
were  the  moments  of  the  unprocessed  speckle  images. 


7.1  Suggested  Extensions  of  this  Work 
The  most  important  theoretical  extensions  to  this  work  would  be 
to  remove  the  assumptions  of  isoplanaticitv  and  stationar ity .  and 
rederive  the  point  spread  function  under  these  new  conditions. 

However,  due  to  the  extreme  difficulty  (or  even  impossibility)  of 
achieving  an  analytic  result  when  these  assumptions  are  removed,  this 
might  well  be  a  point  of  diminishing  returns. 

A  further  theoretical  extension  would  be  to  repeat  the  analysis 
with  the  addition  of  noise;  however,  in  this  case  one  runs  into  the 


problem  of  determining  the  probability  density  function  of  speckles 
plus  noise.  The  speckles  are  lognormally  distributed,  and  the 
lognormal  distribution  does  not  possess  an  analytic  characteristic 
function,  so  the  density  function  of  speckles  plus  noise  would 


necessarily  be  an  approximation. 
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It  is  likely  that  effort  would  be  better  spent  extending  this 
algorithm  to  deal  with  practical  problems  facing  its  use  in  specific 
applications,  since  the  major  theoretical  results  have  been  determined. 
One  such  practical  problem  is  that  of  obtaining  a  good  estimate  of 
*f(u,v)  for  the  Wiener  filter,  and  it  might  well  be  worthwhile  to 
investigate  the  dependence  of  restoration  quality  on  the  accuracy  of 
such  estimates.  Another  problem  is  found  in  astronomy,  where  one  may 
face  extremely  low  light  levels,  say  a  photon  or  two  per  frame. 

Therefore,  in  conclusion,  we  emphasize  that  such  difficulties  are 
likely  to  be  the  most  fruitful  avenues  for  research  on  shif t-and-add 
processing . 
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APPENDIX  A 


THE  LOGNORMAL  DISTRIBUTION 


In  this  appendix,  we  present  some  of  the  basic  properties  of 
the  lognormal  distribution,  as  it  appears  to  be  generally  unfamiliar, 
and  comment  further  on  the  property  that  the  lognormal  probability 
density  function  is  not  uniquely  determined  by  its  moments.  For 
further  information  on  the  lognormal  distribution,  the  reader  is 
referred  to  Aitchison  and  Brown  (1957)  or  Johnson  and  Kotz  (1970). 

The  iuc-paraneter  probability  density  function  of  a  lognormal 
random  variables  X  takes  the  form 


— l—  exp  (-  l‘-°8V  wn 

*  2r  L  “  J 


,  ( x  >  0) 


(  A—  1 ) 


where  log  X  is  normally  distributed,  log  denotes  the  natural  logarithm 

2 

(base  e).  l  is  the  expected  value  or  mean  of  log  X,  and  c~  is  the 
variance  of  log  X.  Following  Johnson  and  Kotz  (1970),  the  nth  moment 
of  X  is 


'  rr vn i  nu  +  l/2n  0 
m  -  E[X  ]  -  e 
n 


(A- 2) 


so  that  the  mean  is 


E(  X ]  *  m^  *  e 


L  +  1/2  0* 
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The  central  moments  may  be  computed  from  the  general  formula 


m  =  - —  ) 

n  (w-l)n/i-  j=0 


( A-  3 ) 


where 


The  distribution  is  unimodal,  with 


mode(X)  =  e" 


;dian(X)  =  e^. 


.Also,  as  o  ->  0,  the  distribution  tends  to  a  normal  distribution. 


Analogs  to  the  Central  Limit  Theorem  also  exist  for  products  of 


random  variables.  Just  as  sums  of  independent  random  variables  are 


approximately  normally  distributed  under  certain  conditions,  so  under 


similar  conditions  are  products  of  independent  random  variables 


approximately  lognormally  distributed.  For  more  specific  statement  of 


the  lognormal  Central  Limit  Theorem  and  for  extensive  results  on 


estimation  of  parameters  for  a  lognormally  distributed  population  or 


data  set,  see  Aitchison  and  Brown  (1957). 


We  now  address  the  property  that  the  lognormal  distribution  is 


not  uniquely  determined  by  its  moments,  a  property  of  interest  in 


Chapter  2.  In  the  usual  case  (Papoulis,  1965),  if  we  know  the  moments 


m  of  a  random  variable,  we  may  determine  its  characteristic  function 


$(u)  by  the  following  expression. 
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where  i  denotes  the  square  root  of  -1.  Inverse  Fourier  transforming 
then  yields  the  probability  density  function.  However,  Hevde  (1963) 
and  later  Barakat  (1976)  have  proved  that  this  result  is  not  unique  for 
the  lognormal  distribution.  Heyde  (1963)  has  presented  another 
distribution  with  density  function 


f(x)  .  — exp  (-  il0Bx;ui 
CxvZfr  2o“ 


}  (1 


+  e  sin 


logx-y)  ]} 


X  >  0 

(A-4) 


where  0  <  £  <  1  and  k  is  a  positive  integer,  which  he  claims  has  the 
same  moments  as  the  lognormal  random  variable  X  with  the  density 
function  p(x)  of  Equation  ( A— 1 ) . 

To  see  this  result,  we  consider  the  difference  in  moments  of 
the  two  distributions 


00 
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I  =  ]  x  p(x)dx  - 
0 


x  f(x)dx  n  *  0,  1,  2,  .  .  . 


-c 

0/2? 


00 

n-1  r  riogx-ul2!  .  r2*k  ,, 

-  x  exp  ^  4 ^ -  >  sin[ — —  (lo 

if  1  l  20"  J  a" 


gx-y) ]  dx. 


We  make  the  change  of  variable 


giving 


-ce  "  |  -l/2v“*ncv  , 2"k  ,  , 
-  l  e  ■  sin  - v  lav 

J  c  '  '  ’ 
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2  2 

-  tuj*1/2ti  a 

-  t.e 


- 1/2(  v-na )"  ,  2mk  ,, 

2  sm[— — yjdy 


because  che  integrand  is  odd.  Therefore  both  distributions  (where 
actually  f(x)  determines  a  whole  family  of  distributions  by  varying  the 
parameters  £  and  k)  have  the  same  moments. 


APPENDIX  8 


evaluation  of  integrals 

In  this  appendix,  we  evaluate  two  integrals  from  Chapter  4 
which  arise  in  the  determination  of  the  overall  point  spread  function 
of  the  shif t-and-add  process. 

We  will  first  evaluate  the  following  integral 

00 

p 

e[i2!i2]  =  j  i2pd2 iii)di2 
o 


here  we  have  collected  terms  according  to  their  dependence  on  io[ 


bus  Equation  (5-1)  becomes 
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2  (l-p“) 


X  exp 


■  -K  (logI2)2  -  2(u2/o2  *  flL-  (logI1-U1))logI2j 


2(1-0") 


( B-2 ) 


e  now  consider  the  integral  portion  only  of  this  expression  and  let 

^  f-[*  (i0§12):  '  2(VC2  +  dogli-u,) ) iogI2j"^ 


2(1-0") 


(3-3) 


c  evaluate  this  integral,  we  make  the  following  change  of  variable  in 
quation  '3-2';: 


t  =  log  I, 


dl„  =  e"dt 
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We  now  evaluate  the  integral  expression  E [ I ^ 1 ^ ] ,  which  is 


needed  in  order  to  determine  the  correlation  coefficient  of  I ^  an 


This  expectation  takes  the  form 
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We  make  the  changes  of  variables 
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First  we  evaluate  the  t  integral 
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by  using  the  same  formula  from  Gradshteyn  and  Ryzhik  ( 1980) 


o  2  Q“  vt" 

exp[-p*x  +  qxjdx  =  exp  — — ,  p  >  0 
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In  this  case,  we  make  the  correspondences 


so  that 
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Substituting  in  Equation  ( B— 6 )  yields 
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which  is  evaluated  as  was  the  t  integral  with 
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